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Preface 


This report is a part of the research and development project on 
aircraft de-icing by the electromagnetic impulse method. This project 
has been sponsored by the Lewis Research Center of the National Aero- 
nautics and Space Administration under grant number NAG-3-284. The 
grant administrator was Mr. John J. Reinmann. For the previous four 
years, many tests had been run in the NASA Icing Research Tunnel and 
three sets of flight tests were performed. These, plus various lab- 
oratory tests, have resulted in a semi -empirical technology for 
designing an electro-impulse de-icing (EIDI) system. 

However, the empirical method is inadequate when a very different 
geometry, material or size is encountered. A computative solution is 
needed which permits prediction of the de-icing effect for a given 
configuration and electrical circuitry. This report, which is princi- 
pally the Ph.D. dissertation of Robert A. Henderson under the direction 
of Prof. Robert L. Schrag, attempts to do the first part of a full 
computer simulation of EIDI. The pressure/time produced by the method 
in this report would be necessary input for a computer code giving 
the structural dynamic response of a given configuration. The con- 
figurations in mind are leading edge portions of aircraft wings, engine 
nacelle inlets and rotor blades. Applications, however, are not 
limited to these applications. 

The author acknowledges the assistance of NASA-Lewis Research 
Center, both for the support of the whole EIDI project at Wichita State 
University and for the opportunity to work at NASA-Lewis during the 
summer of 1985. The assistance of Drs. R. Joseph Shaw, Bill Ford and 
Avram Sidi during that time is gratefully expressed. 
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CHAPTER ONE 


INTRODUCTION 


1. 1 ELECTRO-IMPULSE PHENOMENON AND DE-ICING 

In a paper presented January 8, 1906, at the 24th AIAA Aerospace 
Sciences Meeting in Reno, Nevada, Cl 3 author R. D. Rudich cited air- 
craft icing as the direct cause of four of the thirty two weather 
related fatal accidents reported in the paper. While this may not 
seem like a significant number, the loss of human lives in these four 
accidents could have been prevented if the aircraft involved had been 
able to cope successfully with the icing conditions they encountered. 

The purpose of this dissertation is to present methods of analy- 
zing the electromagnetic aspects of a new method of de-icing both 
private and commercial aircraft. This new de-icing method, referred 
to as Electro Impulse De-Icing (abbreviated EIDI), holds the promise 
of being superior to the present aircraft de-icing methods in terms of 
the energy expended in removing accreted ice [423. 

The EIDI concept is not new. In May 1939, Great Britain issued a 
patent to Mr. Rudolf Goldschmidt covering the basic EIDI mechanism 
C23. However, no commercial development of an EIDI system in the free 
world proceeded from this patent. It has only been vithin the last 5 
years that the unavailability of bleed air from the new high bypass 
ratio engines used on the next generation commercial aircraft has 
caused attention to again be directed to the commercial use of an EIDI 
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system for removing ice from aircraft. 

The simplest practical EIDI system consists of a spirally wound 
coil of rectangular cross section conductor mounted with its axis of 
symmetry perpendicular to the metal surface to be de-iced. An ini- 
tially charged capacitor is discharged through this coil, and the 
resulting magnetic field from the coil's current causes eddy currents 
in the metal. The force exerted on these eddy currents by the c oil's 
magnetic field is initially in such a direction as to cause the coil 
and the metal surface to separate. It is this force that causes ice 
on the metal surface to crack, and subsequently to be removed from the 
surface. 

A considerable body of literature concerned with the electromag- 
netic aspects of a coil placed next to a conducting surface haB accu- 
mulated. Levy [8] and Grover 191 present methods of calculating 
terminal impedances. Dodd and Deeds CIO! discuss both impedance 
calculations and eddy current distributions. Onoe [111 is apparently 
the first researcher to apply a Hankel transformation in the calcula- 
tion of impedances. Onoe's method is further developed and extended 
by El-Markabi and Freeman C41 to include calculation of the force 
between the coil and conductor when the coil current is a sinusoid. 
Bowley et. al. C431 discuss the use of the magnetic vector potential 
in calculating the impulse delivered to the conductor by a transient 
current in the coil. Lewis [441 summarizes the use of the Bowley et. 
al. method for designing a coil to deliver a specific impulse in an 
electro-impulse de-icing installation. 
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1. 2 SCOPE OF DISSERTATION 

An experimental set-up of a prototype EIDI system was assembled 
at The Wichita State University by Dr. Robert Schrag. An experimental 
study of the electro-impulse phenomenon in this system was made by Dr. 
Schrag utilizing field diagnostics methods C3I. Axial and radial 
components of the magnetic field were measured on both sides of the 
rigidly held aluminum "target* plate. The data were then used to 
deduce the total mechanical force versus time, and the mechanical 
impulse strength. Impulse strength was also measured directly with a 
ballistic pendulum. This experimental study provided results which 
will be used to verify theoretical predictions from the mathematical 
model used in this dissertation. 

In this dissertation, the physical phenomena involved in the 
prototype EIDI system of Dr. Schrag 's experiment will be investigated 
analytically and numerically. Specifically, the following tasks will 
be undertaken. 

1. A mathematical model for the total electrical problem will be 
devised. It will employ a transmission line analogy to handle the 
electromagnetic field portion of the system, and a frequency domain 
model for the circuit portion. 

2. The math model will be solved by computer for the specific 
set of conditions that existed in Dr. Schrag 's diagnostics experiment, 
and the calculated and experimental results will be compared. 

1. 3 ORGANIZATION OF DISSERTATION 

Figure 1 summarizes the analysis structure developed in this 
dissertation to predict the behavior of the prototype EIDI system 
briefly described in Section 1. 1. Each of the blocks in this figure 
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represents a stage in the procedure for determining the force-time 
profile on a rigid coil placed next to a fixed conducting plane when a 
capacitor is discharged through the coil. 

In Block 1 of Figure 1, Maxwell's equations are written for the 
coil and metal target of the physical system. A simple model of the 
coil is proposed that stresses the geometrical symmetry of the system, 
eliminating several terms from the Maxwell equations. Chapter 3 dis- 
cusses this modeling of the coil and metal target, and shows that 
application of a Hankel transform to the model equations results in a 
conceptual replacement of the field problem with an infinite set of 
transmission line problems (Block 2). 

In Block 3 of Figure 1, the analysis of the transmission line 
model is performed. The methods of this analysis, for steady state 
sinusoidal conditions, may be found in Chapter 4. Chapters 3 and 4 
provide the basic theoretical developments, with the results of the 
analyses in Chapter 4 applied in Chapter 6 to Dr. Schrag's prototype 
experimental EIDI system, described in Chapter 5. Modeling of the 
circuit used to provide energy to the coil is discussed in Chapter 2 
for the specific experimental system of Chapter 5. 

Once the transmission line model of the coil and target has been 
established in Block 3, Figure 1 shows two possible next steps in the 
analysis. One of these next steps, calculation of the true field 
values (Block 7), requires a knowledge of the coil current. Accord- 
ingly, one proceeds from the transmission line analysis in Block 3 to 
the Block 4 calculation of the coil impedance, required for calcula- 
ting the coil current. Use of the transmission line model for calcu- 
lating the part of the impedance of the coil that is due to the 
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interaction between the coil and target is developed theoretically in 
Section 4.2 of Chapter 4, and applied to the impedance calculation of 
the coil in the prototype EIDI system in Section 6. 2 of Chapter 6. 

Coil impedance obtained in Block 4 is added to the calculated 
skin effect loss resistance in the coil, and with this total coil 
impedance the problem of calculating coil current, using the circuit 
models presented in Chapter 2 , is addressed. This is Block 5, discus- 
sed in Section 6.3 of Chapter 6 for the prototype EIDI system de- 
scribed in Chapter 5. 

Knowledge of the frequency spectrum of the coil current allows 
the calculation of the true field values in Hankel-Fourier space, 
shown in Block 7 and developed theoretically in Section 4.3 of Chapter 
4. Performing an inverse Hankel transformation on these fields 
results in Fourier space (frequency domain) fields. An inverse dis- 
crete Fourier transform applied to the Fourier space fields provides 
the time domain behavior of the electric and magnetic fields, shown in 
Block 8. This procedure is discussed in Section 6. 4 of Chapter 6 for 
the prototype EIDI system. 

Finally, the real space axial and radial magnetic induction 
fields on the coil side face of the metal plate provide the necessary 
information for calculating the total force on the target as a func- 
tion of time, as shown in Block 9. Theoretical development of the 
force equation is in Section 4. 5 of Chapter 4, with Section 6. 5 of 
Chapter 6 describing the numerical implementation of this force calcu- 
lation for the prototype EIDI system. Calculation of the impulse 
delivered to the target is discussed theoretically in Section 4.6 of 
Chapter 4, and its application to the prototype EIDI system is de- 
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scribed in Section 6.6 of Chapter 6. 

Chapter 7 contains a summary of the results obtained, a descrip- 
tion of the original contributions made by the author, and some sug- 
gestions for improving the procedures described in this dissertation 
for modeling the electrical aspects of an EIDI system. 



CHAPTER TWO 


THE ELECTRIC CIRCUIT MODEL 

Modeling of the electronics providing power to the de-icing coil 
in an EIDI system is performed with the desired output from the model 
in mind. Since it is the coil current in conjunction with the 
physical configuration of the coil and wing skin that determines the 
magnetic fields responsible for the forces on the skin, the current in 
the circuit was chosen as the primary variable. 

Testing of a prototype EIDI system began at The Wichita State 
University in 1982. The part of this system that provides power to 
the coil is shown in Figure 5A <page 35). This system was chosen as 
the basic physical configuration for which a circuit model would be 
derived. This circuit model is then analyzed to determine the shape 
of the coil current. A simple circuit model that is suggested by the 
physical system in Figure 5A is shown in Figure 6E (page 47). This is 
the basic circuit model, valid when the clamp diode across the capaci- 
tor is not conducting. 

Justification for modeling the EIDI system as a lumped parameter 
circuit was provided by the experimental verification of the absence 
in the signals present of any significant frequency components having 
wavelengths comparable to the physical dimensions of the system. Be- 
cause of the complex electromagnetic interaction between the coil and 
the skin next to it, it was felt that time domain modeling of the 
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coil's terminal v-i characteristics would be too complex to be of much 
use. Consequently, a frequency domain approach was chosen for 
analyzing the circuit, making the model of the coil a linear frequency 
dependent impedance. 

The presence of the SCR and the clamp diode across the capacitor 
makes the circuit model nonlinear, so that straightforward Fourier 
(frequency domain) techniques are not applicable. This difficulty is 
circumvented by performing a piecewise linear analysis of the circuit. 
In this analysis, the physical circuit is modeled by one of two 
possible circuits, depending on the state of the clamp diode. 

Initially, when the SCR has just been triggered, the diode is assumed 
off (an open circuit) and the model of Figure 6E (page 47) is used for 
analysis. In addition to calculating the current in this circuit, the 
capacitor voltage is also calculated. When this voltage first becomes 
negative, the clamp diode is modeled as coming into immediate forward 
conduction; This diode then acts as a short circuit, resulting in the 
circuit model shown in Figure 6F (page 51). All subsequent (in time) 
circuit calculations are then performed with this model. 

Theoretical justification for modeling both the SCR and the clamp 
diode as simple on-off switches is now provided. However, the 
ultimate justification for such an outright dismissal of the effects 
of both the SCR and the clamp diode in determining the current wave- 
form comes from observing how closely the predicted current waveform 
(using the models that ignore the non-ideal nature of the SCR and 
diode) agrees with the measured waveform. This will be shown in 
Chapter 6. 

In a practical EIDI system, the voltage on the energy storage 
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capacitor prior to its discharge is 1000 to 1500 volts. When the SCR 
is triggered into conduction, its voltage drop is on the order of 1 
volt, which is small compared to the capacitor voltage, and so may be 
ignored in determining the circuit current. In addition to ignoring 
the voltage drop across the SCR, the dynamics of the SCR are also not 
modeled. Such dynamics are, in the EIDI circuit, primarily manifested 
in the failure of the SCR to trigger into instantaneous full forward 
conduction upon initiation of a forward gate current. However, modern 
SCR design techniques C5], C61, 17 ] have resulted in turn-on times 
that are short compared to the time required for a significant change 
in the coil current in the experimental EIDI prototype system. The 
SCR used in the prototype system had a di/dt rating of 800 amps/ 
microsecond, while the maximum observed rate of change in the circuit 
was 15 amps/microsecond. 

A similar consideration of the voltage levels in the circuit re- 
sults in the conclusion that the approximately 1 volt forward drop 
across the clamp diode is not a primary factor in determining circuit 
current. When the diode is off, its transition capacitance is insig- 
nificant compared to the capacitance of the energy storage capacitor. 
Reverse leakage current in the diode is ignored due to the large 
energy storage capacitor in parallel with the diode and the relatively 
small time in which the electrical events of interest take place in 
the circuit. With the diode in forward conduction, the sum of its 
transition and diffusion capacitances are small enough that, to a 
first approximation, they may be ignored in the circuit model. 

In constructing the EIDI prototype experimental configuration, 
care was exercised to minimize parasitics in the circuit. Special low 
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inductance and resistance cable was used to connect the energy storage 
capacitor to the coil. However, both of these cable parameters were 
measured in the prototype system and are taken into account in the 
circuit model. The energy storage capacitor, which was physically 
several units in parallel, used copper strap for wiring connections to 
minimize inductance and resistance parasitics. The ESR (equivalent 
series resistance) and ESL (equivalent series inductance) of the 
capacitors were felt to be small, and are not modeled. If, in a 
particular EIDI installation, these parasitics are not small, they can 
be taken into account by adding lumped elements in series with the 
energy storage capacitor in the circuit model. 

Accurate frequency domain modeling of the coil is the most diffi- 
cult feat in constructing the circuit model, and is discussed in 


Chapter 4. 



CHAPTER THREE 


TRANSMISSION LINE MODEL OF THE FIELD PROBLEM 


3. 1 INTRODUCTION 

The transmission line model of the coil and metal plate is the 
heart of the prototype EIDI system model. It is this model that is 
expected to account for the complex electromagnetic interactions be- 
tween the coil and plate. These interactions help establish the 
impedance presented at the coil terminals, and so are a factor in 
determining the coil current. This current is important in determin- 
ing the force on the plate. 

3-2 GEOMETRY OF A PROTOTYPE EIDI SYSTEM 

Figure 3A shows the profile of the physical coil-metal plate con- 
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figuration we have chosen to model. The coil has the shape of a 
short (h<<R» ) thick walled (R, <<R* ) hollow cylinder with an inside 
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radius R> and an outside radius R* , whose axis is perpendicular to a 
flat metal plate (henceforth called the "target*) of thickness d 
extending to infinity in all radial directions. An air gap exists 
between the coil and target, both of which are assumed rigid and 
stationary in space. 

Most of the coils used in EIDI applications have been round, at 
least prior to a possible bending of the coil to conform to the curved 
leading edge of a wing. The flat geometry of Figure 3A is a reason- 
able model of such a coil, provided that the radius of curvature of 
the coil, after being shaped to conform to the wj.ng, is much greater 
than the coil outer radius R a . 

The assumption that the coil and target remain separated by a 
fixed distance would be reasonable if the initial separation distance 
was much larger than the maximum change in separation distance ob- 
tained when the capacitor is discharged through the coil. Such an 
inequality in separation distances may not hold in a practical EIDI 
installation. Alternatively, the fixed separation distance assumption 
would be justified if the force "impulse" delivered to the plate by 
the coil was so short that the plate acquired only a small velocity, 
with negligible displacement, for the duration of the "impulse". This 
does not generally happen in a practical EIDI system. In fact, a well 
designed installation has the target move from zero to maximum dis- 
placement within the duration of the "impulse*. With ice loading 
present, however, the displacement may be small compared to the ini- 
tial separation distance. 
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3. 3 THE FIELD EQUATIONS 

The approach that appears to be the most fruitful for modeling 
the coil and target in our prototype EIDI system is presented in a 
paper by El-Markabi and Freeman C43. We now describe their approach, 
emphasizing those aspects of the theory that are most appropriate to 
the analysis of the model of the prototype system of Figure 3A. The 
reader is referred to El-Markabi and Freeman C 4 ] for the more general 
theory. 

Current in the coil of Figure 3A is assumed to be entirely phi 

1 

directed. A slight modeling error is introduced by this assumption, 
since a radial component of current actually exists in a practical 
coil due to the spiral winding of the coil. By neglecting this small 
radial current, a model having azimuthal symmetry is obtained. Conse- 
quently, the model shows no phi dependence in any of its field quanti- 
ties, the electric field contains only a phi component, and the mag- 
netic field contains only axial and radial components. 

Because most of the spectral energy of the current in a practical 
EIDI coil is confined to relatively low frequencies, a quasi-static 
situation is assumed. The displacement current term in Ampere's Law 
is ignored. 

With these assumptions, Maxwell's equations written for Figure 3A 




Partial differentiation with respect to time can be eliminated from 
these equations by taking the Fourier transform. Differentiation with 
respect to time is then replaced with multiplication by the opera- 

tor. Performing this transformation, we obtain the following equa- 
tions (note that we have not introduced new symbols for the trans- 
formed field quantities, but have instead simply replaced their time 
variable argument t with the Fourier transform variable omega) 


3 * dr- * ’ > ' 
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Equation (1) contains only the phi component of E and the radial 
component of H. By combining equations (2) and (3) in such a fashion 
as to eliminate the axial component of H, we would have another equa- 
tion containing only the phi component of E and the radial component 
of H. This elimination yields 
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At this point, we introduce a mathematical tool that is of con- 
siderable use here. This tool is the Hankel transform of order n, 
defined by C 1 2 3 

S j rT^Xr^lr ^ F {\) 

J a 


with the inverse transform 
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Although the introduction of an additional transform introduces 
difficulties of its own, it provides an even greater amount of simpli- 
fication, and makes possible the transmission line model of the coil- 
target electromagnetic field problem. 

If a Hankel transform of order 1 is applied with respect to the 
variable r in equation (4), the following result is obtained till, 

C123 (note that we have again not introduced new symbols for the 
transformed quantities, simply replacing their spatial variable argu- 
ment r with the Hankel transform variable lambda) 

9l4 r 








The advantage of equation (5) over equation (4) is that equation (5) 
contains partial derivatives with respect to only the z coordinate, in 
effect becoming an ordinary differential equation (if one allows ’con- 
stants of integration* for such an equation to be arbitrary functions 

l 

of all variables except z). If the same Hankel transform is applied 
to equation (1), we obtain 


at 




((o') 


Equations (5) and (6) are two coupled ordinary differential equations 
that are recognizable as the canonical transmission line equations. 
The solutions of these equations are well known C 1 3 1 , (14), C15J. 


3. 4 DEVELOPMENT OF THE TRANSMISSION LINE MODEL 

Because equations (5) and (6) are identical in form to the equa- 
tions describing voltage and current on an ordinary transmission line, 

A A. 

we now introduce the symbols V and I (each of which is a function of 
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position z along the coil axis, the Fourier variable omega, and the 
Hankel variable lambda) to represent the phi component of E and the 
radial component of H respectively. Then equations (5) and (6) become 
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Differentiating 

equation (7) with respect to z, ve have 
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equation (8) into this last equation yields 
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Similarly, we differentiate equation (8) to obtain 

a a v . ax 
a a* = 6*^* a* 


Substitution of equation (7) into this equation yields 



Nov define the complex propagation constant gamma in Fourier-Hankel 
space as 


X = -V X* 4- <3- 

so that the general solutions to equations (9) and (10), assuming com- 


plex sinusoidal time variation of V and I, may be written as 

-*a- 
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A. -a. A A 

where !<>♦, Io- and V 0 ., V 0 - are complex phasor functions of the Hankel 


variable lambda (we have suppressed the complex sinusoidal time varia- 
tion ev in these solutions). 


Mote that we are conceptually dealing with an uncountably infi- 
nite set of transmission lines. Each transmission line is associated 


with a different value of lambda, where lambda is non-negative. The 
physical origin of this infinite number of transmission lines lies in 
our "compressing" the infinite radial variation in the real space 
field quantities H, and s at a given axial coordinate z into vari- 
ables V and I localized to a single "point" (the corresponding z 
coordinate) on a Hankel space transmission line. It then takes an 
infinite number of transmission lines to account for the infinite 


number of possible values of the radius in the real space problem. 

To derive an expression for the characteristic impedance of one 
of these Hankel space transmission lines, consider the case of a line 
having only a single frequency complex sinusoid traveling in the 
direction of increasing z. In agreement with ordinary transmission 
line theory, the negative sign in the exponents in equations (11) and 
(12) denotes propagation in the direction of increasing z. Then 
equations (11) and (12) become 
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Substitute the expressions for the current and voltage from equations 
(13) and (14) into equation (7) 
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From this equation, we obtain the defining expression for the charac- 
teristic impedance of a Hankel space transmission line as 
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This result differs by a negative sign from the expression given 
in El-Markabi and Freeman C4] for the line's characteristic impedance. 
One could argue that the negative branch of the square root function 
in the expression for the propagation coefficient gamma should be 
chosen, which would eliminate the negative sign in the characteristic 
impedance. However, calculations performed with a negative character- 
istic impedance results in predictions that are essentially in agree- 
ment with experimental measurements, whereas the use of a positive 
characteristic impedance predicts results that are not in agreement 
with experiment. Alternatively, one could propose that it is the 
positive sign that must be chosen for the exponents in equations (11) 
and (12) to correspond to propagation in the direction of increasing 
z. Such an assumption also results in predictions that are not in 
agreement with experiment. 

Transformation of real space current sources into Hankel space is 
discussed in El-Harkabi and Freeman C4). They show that a disc of 
azimuthally directed uniform surface current of phasor value I, with 
an inner radius Ri and an outer radius R», located on the z-axis, 
becomes a sinusoidal current source with phasor value 
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connected in parallel with the Hankel apace transmission line. The z 
coordinate of this current source is the same as the z coordinate of 
the real space surface current from which it arose. 

We have chosen three of these ideal azimuthally directed discs of 
uniform surface current to model the current distribution in the coil. 
The inner (outer) radius of the disc is equal to the inner (outer) 
radius of the coil. The middle disc is located at the axial center of 
the coil, while the two outer discs are each spaced out a distance of 
one-third of the coil's width from the center of the coil. See Figure 
3B. Clearly, some error occurs here due to the localization of the 
model's current to these discs. 



FIGURE 3B 

CURRENT DISC MODEL 
OF COIL 


Figure 3C shows the Hankel space transmission line configuration 
corresponding to the three current sheet model of the coil next to the 
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target. All of our Hankel space calculations are based upon this 
model. 
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FIGURE 3C 

HANKEL SPACE TRANSMISSION LINE MODEL 
OF COIL AND TARGET 



CHAPTER FOUR 


ANALYSIS STRUCTURE 


4. 1 INTRODUCTION 

In this chapter, the transmission line model of the coil and 
target developed in Chapter 3 and shown in Figure 3C is examined 
mathematically using conventional transmission line theory to predict 
the Fourier-Hankel space behavior of the model when it is in the 
sinusoidal steady state condition. The analyses performed in this 
chapter will assume each of the three current sources in Figure 3C has 
unit phasor value. In view of the analogy developed in Chapter 3, we 
will freely use the terms current and voltage in connection with the 
transmission line model to stand for the Fourier-Hankel transforms of 
the radial magnetic induction B, (z, r, t) and the azimuthal electric 
intensity E^z, r,t) respectively. Throughout this chapter, reference 
should be made to Figure 3C for insight into the equations written and 
for identification of the variables used. 

4.2 METHOD OF CALCULATING COIL IMPEDANCE 

Although the part of the coil impedance that is independent of 
.the properties of the ribbon conductor used to wind the coil is calcu- 
lated in Hankel space, as shown in this section, one must begin the 
derivation of the expresion for the impedance in Fourier space. Only 
in the frequency domain does the concept of impedance make sense. We 
will calculate the coil impedance at a frequency omega by exciting the 
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coil with a current source of phasor value 1, and finding the phasor 
voltage response at the coil terminals. The coil impedance is then 
numerically equal to this voltage response. 

Qur real (Fourier) space model of the coil is three ideal discs 
of azimuthally directed uniform surface current. We need to relate 
the electric fields induced in the different parts of this model to 
the single voltage that exists between the coil terminals. We approx- 
imate the voltage appearing at the coil terminals by the average of 
the voltages induced at the locations of the three discs of the model. 
Symbolically, we have 


? = 


V_ 

A/ 

X. 


\ 




'V 

where V* represents the voltage induced at the location of the disc 
with axial coordinate z* . 

Because of the symmetry of our model, the voltage induced on a 
circular path (centered on the z axis) of radius r and axial coordi- 
nate z* will be given by 


, r } ui) 

The average voltage induced on the infinite number of circular paths 
between r=Ri and r=R* is 



Multiplication of this expression by N, the number of turns in the 
actual coil, yields 



r 


r ; 



2k 


which is the model -predicted voltage induced on an infinitesimally 
thin N turn coil having an inner radius Ri and an outer radius R 2 
located at z k . Averaging these induced voltages over the three discs 
results in the phasor coil voltage 


V 


* J 1 r Bo & 4>^0 *) M 


V - z 


Because this expression is numerically equal to the coil impedance, we 

A/ 

will replace the symbol V by Z in further appearances of this expres- 
sion. 

The E fields that appear in (1) are frequency domain (Fourier 
space) E fields. In terms of the Fourier-Hankel space E fields, we 
can write these Fourier space E fields as 


OO 

© 

vhere we have again used the list of arguments to distinguish between 
different functions. Specifically, E^z* , r, ) is the Fourier trans- 
form of the real space azimuthal electric intensity E^,(z k ,r, t), and 
E^z* , X , o> ) is the Fourier-Hankel transform of the real space azimu- 
thal electric intensity E^Zn , r, t). Substituting (2) into (1) yields 
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Interchanging the order of integration. 
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where ^ 

^ 1 ^ T ' (x ^ r Jr 

is the Hankel transformed current corresponding to an infinitesimally 
thin H turn coil of inner radius R, and outer radius R* , carrying a 
phasor current of strength 1/3. 

Expression (3) contains three Fourier-Hankel space electric 
fields (transmission line voltages) that conceptually arose from Four- 
ier space phasor currents of value 1/3 on each of the three discs in 
Figure 3B. The "scaling factor* that relates Fourier-Hankel space 
fields due to unity Fourier space phasor coil current to Fourier- 
Hankel space voltages and currents due to unity phasor current in each 
of the three current sources in Figure 3C is K'l X). Then in terms of 
the Fourier-Hankel space voltage E^z, , X , <») due to unity phasor 
current in each of the Figure 3C current sources, (3) becomes 

2.-- 4*5 K ^ 

o 

Expression (4) is the result that ve use to calculate the terminal 
impedance of the coil. 

4.3 RADIAL MAGNETIC INDUCTION AND AZIMUTHAL ELECTRIC INTENSITY 
4.3.1 Target Surface Facing Coil 

We begin by calculating the characteristic impedance and the com- 
plex propagation constant of the air transmission line and the metal 
transmission line, referred to hereafter as the air line and the metal 
line respectively, using the results derived in Chapter 3. For the 
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air line, the conductivity <r* is equal to zero, and so 




X titJiur 


<T/' 


- ^ 


- V a a ~ = ^ 

For the metal line, no simplification is possible and we have 

— 




+ i iJi> C"* 


r <3“/ 


According to the analogy developed in Chapter 3, calculation of 
the total current and voltage at z=0 in Figure 3C is equivalent to 
calculation of the radial magnetic induction B, and the azimuthal 
electric intensity E^ on the coil side face of the target. We perform 
this current and voltage calculation by replacing the transmission 
line configuration for all z>0 with its equivalent impedance. Note 
that both sides of the metal line are connected to infinite lengths of 
air transmission line with characteristic impedance Z. . The equiva- 
lent impedance looking to the right at z-0 in Figure 3 C is 


2 1 °) 


(*„j) 


' TA 
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Using this equivalent impedance, the current reflection coef- 
ficient at z=0 as seen from the air line is 

' l a + Z(o) 

and the total current at z=0 is then 

A. A A A / \ 

I| O I « I = line * Ir • f I • et id = 

A 

in terms of the incident current I,„ 0 , which must be due solely to the 
three current sources. Expression (6> may be simplified as follows. 
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Substituting (5) into (7) yields 
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which, upon replacing Z« and Z. by their defining expressions gives 
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Now consider the current source in Figure 3C closest to the tar- 
get. This source produces an incident current on the metal line at 
z=0 given by 

3 . 


A 


. A."** 
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where the factor 1/2 comes from the equal division of the unity phasor 
amplitude current to both sides of the air line containing the current 
source. Total current incident at z=0 is by superposition given by 


A 

C. 

i - v a * i j u 

A. 
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Substituting (9) into (8) yields 
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for the total steady state current I«„«.i at z = 0 due to unit strength 
complex sinusoidal excitation in each of the three current sources. 
Multiplying (10) by K y < X )I («•> ), where 1(^0) is the Fourier transform 
of the current in the coil (and multiplication by K '(. X ) transforms a 
unity Fourier space current into the corresponding Fourier-Hankel 
space current), yields the desired result 
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for the Fourier-Hankel transform of the radial magnetic induction B, 
on the coil side face of the metal target when the coil current 
spectrum is I(oO). 

Calculation of the voltage is almost identical to the calculation 
of the current performed above. The voltage reflection coefficient 
is the negative of the current reflection coefficient, 

Z(o) - 

fs ~ ^ (a ) 4- * 


Total voltage at z=0 is 


A 

V. 
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in terms of the incident voltage V, 


Expression (12) may be simpli- 


fied as 
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Substituting (5) into (13), 
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which, upon replacing Z, and Z. by their defining expressions, gives 
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Total incident voltage V t , e at z=0 is given by 

-\<V , -\Wa 
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using (9) for I iI>e . Substituting this expression for V,, e into (14) 


yields 


V/ 


+c4al 


\<JU /. . -Xlo/3 -3^\n/z\ 

<2. 0 ( I + e -i-e.. J 

Ite'nii ( 


(15) 
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for the total voltage at z=0 due to unit strength complex sinusoidal 
excitation in each of the three current sources. Multiplying (15) by 
K / (X)I(oo) yields the desired result 
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for the Fourier-Hankel transform of the azimuthal electric intensity E^ 
on the coil side face of the target when the coil current spectrum is 
1 ( 60 ). 

4. 3. 2 Target Surface Opposite Coil 

Knowing the total voltage and current at z= 0 (derived in Section 
4.3.1 and given by (16) and (11) respectively) allows a simple trans- 
mission line inverse chain matrix calculation of the total voltage and 
current at z=d, which are the Fourier-Hankel transforms of the azimu- 
thal electric intensity E^ and the radial magnetic induction B r on 

the surface of the metal target opposite the coil. Denoting these 

>\ 

transforms of the fields on the coil side face of the target by V. 

A 

(given by (16) above) and I n (given by (11) above), we have 

Yp. ~ Yi c ^ /l ~ (17) 

I, * -V„ („) 

WV| 

^ A. 

where V, and I, denote total voltage (electric intensity) and current 
(magnetic induction) on the surface of the target opposite the coil. 
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4. 4 CALCULATION OF THE AXIAL MAGNETIC FIELD USING THE 
TRANSMISSION LINE MODEL 


The procedure that we use for calculating the force on the target 
is described in the next section, and requires both the radial and the 
axial components of the magnetic induction on the target surface next 
to the coil. The radial component of the magnetic intensity is avail- 
able from the transmission line model by performing inverse Hankel and 
Fourier transformations on the calculated transmission line current at 
z=0, as discussed in Section 4.3.1. By using equation (2) in Chapter 
3, it is also possible to calculate the axial component of the magnet- 
ic intensity (which does not have a transmission line analog) from the 
transmission line voltage, as shown below. 

The azimuthal electric field is given in Fourier space by 

Substituting (19) into equation (2) from Chapter 3, we have 




- I 


trr 







Interchanging the order of partial differentiation and integration in 


( 20 ), 








Using the identity 


(21) becomes 
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Equation (22) states that the axial component of the magnetic intensi- 
ty H, in Fourier space is given by the zero order inverse Hankel 
transform of the transmission line voltage at the corresponding z 
coordinate multiplied by lambda, with the result divided by j 
Inverse Fourier transformation of the right hand side of (22) yields 
the desired time domain axial magnetic intensity. 

4. 5 FORCE BETWEEN TARGET AND COIL 



Our procedure for calculating the total force between the target 
and the coil utilizes the Maxwell stress tensor, and is performed in 
real space (using the time domain fields). Stratton C41J shows that 
the total force F transmitted by a time varying electromagnetic field 
across a closed surface S is given by 


F = ^ [ £ (£»ft\t§L + jr(e.*n)& - yr) ^ 
s 

where ^ is a unit vector normal to the surface. We take as our closed 
surface the plane z=0 (the coil side face of the target, "closed" at 
infinity). Since n=z“, this integral reduces to 

^ = I 5 + + 

r=o 4> = o ' 


The total force tending to separate coil and target is just the z 
component of this force, 
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F^ = J " &r ) r - 'ir &4> r ^ r ( A ^) 

Our calculations of the fields E ^ , B, , and B, in the prototype 
experimental EIDI configuration described in Chapter 5 showed that the 
second integral in (23) is totally insignificant compared to the first 
integral. Some feeling for why this is true can be obtained by com- 


paring the constants that multiply each of the integrals in (23). 
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Accordingly, we approximate (23) by 

- f j ^ (**) 

o 

Equation (24) is the result that we use to calculate force versus 
time. 


4.6 IMPULSE DELIVERED TO TARGET 

The impulse delivered to the stationary target is by definition 
given by the integral 

,«) 

© 

where F, (t) is calculated using (24). 



CHAPTER FIVE 


A SPECIFIC SYSTEM EXAMPLE, INCLUDING EXPERIMENTAL RESULTS 


5. 1 DEFINITION QF THE SYSTEM 

The prototype EIDI system constructed at The Wichita State Uni- 
versity, and the results of the tests made on that system, have been 
described in detail by Dr. Robert Schrag C 3 3 . Most of the material in 
this Chapter has been excerpted from Dr. Schrag \s paper. 

Figure 5A shovs the prototype EIDI energy discharge system, omit- 
ting the capacitor charging circuit and the thyristor firing circuit. 
Two identical pulsing coils were operated in series, because that was 
the arrangement used in most of the de-icing tests. However, only one 
of the two coils was utilized in the coil-target assembly, which is 
detailed in Figure 5B. 

The effective gap between the coil (copper) surface and the near 
surface of the target was .078 inch. A .032" thick 2024 T3 Aluminum 
disc was used as the target. Diameter of this disc was 5 inches. Two 
.05" thick phenolic spacer plates were used, one to maintain a fixed 
distance between the coil and the target, and the other to maintain 
the distance between the target and the rigid wooden support that 
prevented motion of the target. These plates could be removed, and a 
special magnetic field measuring plate (described in Section 5.4) 
inserted in their place in order to make measurements of the magnetic 
induction close to the surface of the target. Each coil consisted of 
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30 turns of .024" X .188" rectangular copper wire spirally wound in a 
single layer from an inner radius of .125" to an outer radius of 1". 
The initial capacitor voltage utilized for the experimental study was 
400 volts. 

5. 2 COIL IMPEDANCE MEASUREMENTS 

Impedance measurements on the coil next to the target were made 
using an impedance bridge. The inductance determined by these mea- 
surements, for several frequencies, appears in the table below. Impe- 
dance measurements were also made on the coil when the metal target 
was removed. The real part of the impedance measured on the coil 
without the target in place was subtracted from the real part of the 
impedance measured on the coil with the target in place. This in- 
crease in resistance due to the target is given in the table below. 

RESISTANCE INCREASE AND INDUCTANCE - COIL AND HETAL TARGET 

Frequency Inductance Resistance 

( Hertz) ( Microhenries ) Increase (Milliohms) 


500 

18.6 

8.2 

1000 

17.0 

23.0 

2000 

14. 6 

48.0 

4000 

12. 6 

77.0 


5. 3 CURRENT WAVEFORM 

Initially, current in the coil was measured indirectly by measur- 
ing the voltage across the .001 ohm non-inductive resistor in Figure 
5A. Difficulties with this approach prompted the purchase of a cur- 
rent transformer from Pearson Electronics, Inc. Using this transfor- 
mer and a storage oscilloscope, we observed the current shown in 




3? 


Figure 5C. 



FIGURE 5C 
COIL CURRENT 

5.4 MAGNETIC FIELD MEASUREMENTS 

A magnetic field measuring plate was constructed in the manner 
illustrated in Figure 5D. Shallow concentric grooves were cut into 
both sides of a .05 inch phenolic disc, with radius increments of .2", 


TYPICAL GROOVES 
WITH SINGLE 
TURN! WIRE 
LOOPS 


CHANNELS WITH 
TWISTED LEADS 



7 DIAMETER 
PHENOLIC DISC, 
. 05 " THICK 


SOLDER TAES 


FIGURE 5D 

PARTIAL ILLUSTRATION 
FIELD MEASURING PLATE 
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starting at r = .2* and ending at r=2.Q". Single turn loops of .006* 
diameter wire were then cemented into these grooves, and their twisted 
leads brought out to solder tabs through radial channels. 

For measuring the fields on either side of the target, the mea- 
suring plate simply substituted for the corresponding phenolic spacer 
plate in Figure 5B. A measurement of the axial flux density was 
derived from the induced voltage in any two neighboring loops connec- 
ted in series opposition. For the two loops illustrated in Figure 5D, 


for example. 




j Ax. 

* <*■/ - r ,^ 


where B, is in teslas, r is in inches, and V is in volts. This value 
is the average axial flux density over the area between the two induc- 
tion loops. In the further use of this result, we will assume the 
flux density to apply at a radius midway between the two loops. 

To measure the radial component of flux density at any radius, 
the front and back loops at that radius are connected in series oppo- 
sition, and calculations are made from 

o 

2-rr K 


%(*) = 1*50 


where r is the radius of the two induction loops and h is their 
separation, both in inches. 

Plots of the magnetic induction fields obtained from these tests 
with the target in place are shown at selected radii in Figures 5E, 

5F, 5G, and 5H. All B, data showed an anomolous behavior (irregu- 
larities) at r=.4" relative to r=.6". A separate check was made, in 
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which the plate was reversed (interchanging the two sides). This 
produced a reversal of the irregularities, so the effect was probably 
due to an inaccuracy in the construction of the plate. 

5. 5 MEASUREMENT OF IMPULSE TO TARGET 

The final step in the experiment was the measurement of the 
impulse delivered to the target when the capacitor was discharged 
through the coil. This was done indirectly, using a ballistic pendu- 
lum. The impulse so measured was approximately . 012 lb-sec. 




CHAPTER SIX 


COMPUTER ANALYSIS ON THE SYSTEM EXAMPLE 


6. 1 INTRODUCTION 

Theoretical background for the numerical computations described 

in this chapter were developed in Chapters 3 and 4. This chapter will 

r 

concentrate on a description of the numerical methods used to imple- 
ment the theoretical development contained in these earlier chapters 
to predict the performance of the prototype EIDI system described in 
Chapter 5. The sections that follow are arranged in the order of the 
Figure 1 Analysis Flow Diagram, beginning with Block 4 of that dia- 
gram. This is the order in which the computations were actually 
accomplished. 

All computations were performed in FORTRAN IV on an IBM 370. 
Source code was compiled with the IBM furnished G level compiler. 

6.2 COIL IMPEDANCE 

Equation (4) of Chapter 4, repeated as equation (1) below, is the 
coil impedance predicted by the transmission line model. Numerical 
evaluation of the integral in (1) is discussed in this section. 

o ~ 

Two common methods are in use to allow quadrature of an improper 
integral such as the Integral above £21]. In the first method, a 
transformation of variables is performed prior to construction of an 
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algorithm for estimating the value of the integral. This transforma- 
tion is chosen in such a manner that the new integral has finite 
limits. The second method simply replaces the infinite upper limit 
with a finite upper limit, selected such that the part of the integral 
thus ignored contributes little to the true value of the integral. 

Such a method can be employed only if the integrand decays sufficient- 
ly rapidly as the variable of integration increases. Since the con- 
trolling factor in the leading behavior of the asymptotic expansion, 

— V q 

as lambda tends to infinity, of the integrand in (1) is e 3 , with g 
a constant, the second method was chosen for use in the numerical 
evaluation of (1). 

A large amount of high quality mathematical software is available 
today [221, [23]. Because of the complexity and cost of writing 
quality mathematical algorithms, most numerical analysts suggest that 
complex scientific calculations be performed using algorithms written 
by experts [24], [25]. An B panel adaptive Newton-Cotes algorithm 
entitled QUANC8, described in [25], was chosen to perform the quadra- 
ture in (1). A user of QUANC8 may select the relative and absolute 
error performances desired, and the program then attempts to estimate 
the value of the integral within the selected error criteria. One of 
the parameters in the subroutine QUANC8 is an output variable that 
contains an estimated error bound on the returned value. 

Direct computer evaluation of (1) consumes a large amount of CPU 
time, and is consequently expensive. This is due mostly to the ap- 
pearance of the factor ( K ' ( X)]* in the integrand. (Note that K / (X) 
is defined by an integral containing a Bessel function. This makes K' 
oscillatory, shown in Figure 6A. > It was possible to decrease this 



44 


cost considerably by calculating !<'( X) using QUANC8 and fitting a 
cubic spline function to the calculated K'< \> for use in evaluating 
(1). The algorithms used for generating the spline function coeffic- 


FIGURE SA 
OSCILLATORY 
NATURE OF K '< X > 


ients and for evaluating the spline function at a given argument are 
entitled SPLINE and SEVAL, respectively, and are described in [25]. 

Once a spline function approximation for K'(X) was available, 
QUANC8 was used to individually estimate the real part (the coil 
resistance increase due to the metal target) and the imaginary part 
(the coil reactance) of (1) at selected frequencies. Appendix A 
contains a listing of the complete program, with the subroutines 
referred to, for evaluating the real part of (1). A simple modifica- 
tion of this program allows the imaginary part of the impedance to be 
evaluated. 

Initial estimates of the real and imaginary parts of the integral 
in (1) were obtained using an upper limit of 680. Following this, the 
programs were run again with an upper limit of 240. Little difference 
was seen in the results of the two calculations, leading to the con- 
clusion that the coil resistance and inductance estimates were good. 






Figure 6B shows the effects of the upper limit of integration in 
evaluating the imaginary part of the integral in (1) for three differ- 
ent frequencies. 



Figures 6C and 6D on the following page show the resistance in- 
crease and inductance calculated as described above. Measured values 
of these parameters are also shown for comparison. 

6.3 CURRENT WAVEFORM 
6. 3. 1 Introduction 

As discussed in Chapter 2, three factors preclude a simple calcu- 
lation of the coil current in the EIDI prototype system. These three 
factors are the nonlinear diode and SCR, and the presence of the metal 
target next to the coil. 

The effects of the target on the coil were taken into account by 
calculating the coil resistance increase and inductance at several 
frequencies, as described in the previous section. These two frequen- 
cy dependent parameters were then approximated with cubic spline 
functions, again using the subroutine SPLINE. This provided the 
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ability to calculate computationally inexpensive yet sufficiently 
accurate coil impedances for use in frequency domain circuit calcula- 
tions. 

6. 3. 2 Current Before Clamp Diode Conduction 

Although the circuit is nonlinear due to the diode and SCR, it is 
amenable to treatment in a piecewise linear fashion. The frequency 
domain model of the circuit, valid between the time the SCR is init- 
ially triggered and the time the capacitor voltage becomes zero, is 
shown in Figure 6E. In this Figure, the EIDI system capacitor is 
modeled as an ideal discharged capacitor in series with a voltage 



FIGURE 6E 
FREQUENCY DOMAIN 
CIRCUIT MODEL 
CLAMP DIODE OFF 


source. This voltage source has a value of zero volts for all time 
prior to t=0. At t=0, it instantly rises to the amplitude and polar- 
ity V 0 of the initial voltage on the actual capacitor, modeling the 
triggering of the SCR into instantaneous full forward conduction. 
Voltages and currents calculated from this model are good approxima- 
tions to their corresponding quantities in the physical EIDI prototype 
circuit as long as the diode in parallel with the capacitor is not 


conducting. 
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The frequency dependent resistor appearing in Figure 6E is a com- 
posite lumped model of several loss mechanisms in the physical cir- 
cuit. It includes the loss in the coil due to the presence of the 
metal target next to the coil (calculated as described in the previous 
section), resistance of the ribbon conductor used to wind the coil 
(corrected for skin effect), and the resistance of the cable used to 
connect the coil and capacitor. Resistance in the cable was modeled 
as frequency independent, with a measured D. C. value of .054 ohms. 

Coil winding resistance was calculated from 

where 

Roc = D. C. coil resistance = .0235 ohms 
h - coil thickness = .00477 meters 
$ = copper skin depth = l uq- 

f = frequency in Hertz ' 

£ = copper permeability = 4^ X 10* 7 henries/meter 
<r* = copper conductivity = 3. 48 X 10 7 mhos/meter 

Resistance R #c applies to both the coil next to the target and the 

idler coil (see Figure 5A). 

Similarly, the frequency dependent inductance L in Figure 6E 
arises from more than one source. First, there is the inductance of 
the coil next to the target (calculated in the section above). Sec- 
ond, the idler coil and cable connecting the capacitor to the coil had 
a combined measured inductance of 23 microhenries, which was assumed 
to be independent of frequency. As these two inductances are in 
series in the circuit, they are added together to obtain a value for 
L. 


The leakage reactance and measuring circuit impedance reflected 
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into the actual circuit by the current transformer used to measure the 
current were felt to be too small to be significant, and were not 
included in the Figure 6E model. 

By inspection, the Fourier transform of the current in Figure 6E 
is given by 
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Mote that the delta distribution multiplies its own argument, and 
consequently contributes nothing to the time domain current i(t). 
Equation (2) is the basis for the current calculation, valid until the 
capacitor voltage becomes zero and the clamp diode in parallel with 
the capacitor begins conducting. For simplicity and cost considera- 
tions, an inverse discrete Fourier transform (IDFT) (361 was chosen 
to approximately invert I ( <o ) . To sufficiently minimize the effects 
of the IDFT approximation to the continuous inverse Fourier transform, 
a 1024 point transform with a sample time of 5 microseconds was cho- 
sen. This makes the folding frequency 100kHz, which is considerably 
above any significant frequencies measured in the spectrum of the 
current in the prototype EIDI system. The 200kHz frequency window of 
the IDFT is sufficiently large that "windowing* effect errors in the 
time domain response are not too great. (These errors are primarily 
manifested as Gibbs' ripples on the intitial current calculated in 
Section 6. 3. 3, described in Section 6. 3. 4 and shown in Figure 6G. ) 
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Appendix B contains a listing of the program used to calculate 
the coil current in the manner described above. 

In order to determine how long a time this calculated current 
approximates the actual circuit current, a second calculation was per- 
formed using the circuit model in Figure 6E. The model predicts a 
frequency domain voltage corresponding to the physical EIDI capacitor 
voltage of 
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The term involving the delta distribution again contributes nothing to 
the inverse Fourier integral yielding the time domain voltage, and is 
dropped from consideration to yield (3). A 1024 point 5 microsecond 
sample time IDFT was used to approximately calculate the inverse 
Fourier transform of (3). Output from this program shoved that the 
capacitor voltage would slew negatively through zero volts at t=279 
microseconds. This is the time at which the diode across the capaci- 
tor is modeled as coming into conduction and acting as a short cir- 
cuit. Beyond this time the model of Figure 6E is no longer valid, and 
consequently the current predicted by the model is no longer correct. 


6.3.3 Current After Clamp Diode Conduction 

Once the diode comes into conduction, the model of the EIDI 


circuit for all future time is as shown in Figure 6F. The inductance 
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FIGURE 6F 
FREQUENCY DOMAIN 
CIRCUIT MODEL 
CLAMP DIODE ON 


and resistance in this circuit have the same physical significance 
that they had in Figure 6E, and their values are calculated for any 
desired frequency using the same algorithms. For computational con- 
venience, time is * reset* to zero in this circuit, even though the 
circuit does not come into existance until t=279 microseconds in the 
Figure 6E circuit. Since the current in the coil is not initially 
zero, the circuit's magnetic energy storage is modeled by including a 
DC current source in parallel vith a frequency dependent inductance. 
This current source is zero for all time prior to t=0, and for all 
future time has a D. C. value I„ equal to the value of the current in 
the Figure 6E circuit at the time the capacitor voltage became zero. 

Current in the Figure 6F circuit is given in the frequency domain 

by , \ 

\lj Lfcu ) 

X(o>) = +. 

\ U 

* X * lei L(<A + K/gj} 

U 

Note that the term containing the delta distribution once again con- 
tributes nothing to the inverse Fourier transform of I (to), and is 


dropped in (4). An IDFT was used to calculate the approximate inverse 
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Fourier transform of 1(^0) given in (4), in exactly the same manner 
that the expression for I(<^o) in <3) was inverted. The program is 
very similar to the one used in the Section 6.3.2 current calculation, 
contained in Appendix B. 

6. 3. 4 Combining Pre- and Post Clamp Diode Conduction Current 

Because of the presence of the current source in parallel vith 
the inductor in Figure 6F, the current in this circuit model is dis- 
continuous at t=0. Consequently, the current cannot be bandlimited. 

An inherent assumption in the use of the IDFT to approximately calcu- 
late a continuous inverse Fourier transform is that the signal is 
bandlimited. The time domain current calculated by the IDFT showed a 
small amount of ripple during the first 100 microseconds due to the 
current spectrum not being bandlimited. Prior to "joining* in time 
the current predicted by the Figure 6E model with the current predict- 
ed by the Figure 6F model, this ripple was eliminated by graphically 
choosing current values lying close to the IDFT predicted values that 
joined smoothly with the current from the Figure 6E model. Figure 6G 
illustrates this procedure. Current ripple amplitude was 54 amps peak 
to peak at t=10 microseconds, decaying to 1 amp peak to peak at t=75 
microseconds, with time measured in the Figure 6F model. Figure 6H 
shows how good the agreement is between this piecevise linear model 
predicted current and the current measured in the prototype EIDI 
circuit. 

With the model predicted current now available at 5 microsecond 
intervals from t=0 to t=5. 12 milliseconds, a 1024 point 5 microsecond 
sample time DFT was performed on the current samples to estimate the 
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Fourier transform I (oo) of the model current. This transform is 
needed for use in calculating the magnetic induction fields from the 
Hankel space transmission line coil-target model, as discussed in the 


next section. 
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6. 4 MAGNETIC INDUCTION 


6. 4. 1 Introduction 

Knovledge of the radial and axial magnetic induction on the coil 
side face of the target is required for calculation of the total force 
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on the target. Although conceptually simple, the numerical calcula- 
tion of these two fields presented more computational difficulties 
than were encountered in the combined total of all other calculations 
performed. Initial attempts at calculating the magnetic induction, 
performed using single precision arithmetic and using professionally 
written quadrature routines, produced results that were incorrect by 
orders of magnitude. 


6.4.2 Radial Magnetic Induction on the Coil Side Face of the Target 

Expression (5) shows the iterated improper integral that is to be 
evaluated numerically to calculate the time domain radial magnetic 
induction B, <0, r, t>. This integral is the inverse Fourier-Hankel 
transform of the total current at z=0 given by (11) of Chapter 4. 






os - — 


oo 




[i +- X Jx 


l + 


X 








-4 


4- 


X - o 


r/' 


*3/J» ( i') 


Quadrature of iterated integrals is almost always difficult C 21 ] . 
One of the most common methods of evaluating such integrals, the Monte 
Carlo method, could not even be considered for use in (5) due to the 
enormity of the very expensive and random complex function evaluations 
required by such an approach. Furthermore, the integrand in <5) is 
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extremely oscillatory, having both positive and negative complex val- 
ues, due to the Bessel function in the inverse Hankel transform, the 
function K '( X), and the complex exponential in the inverse Fourier 
transform. Quadrature of such integrands is all but impossible using 
Monte Carlo methods due to extreme smearing C25], [26]. It is the 
oscillatory nature of the integrand in (5) that makes its evaluation 
difficult. 

After much trial and error, a workable approach to quadrature of 
the integral in (1) was obtained, and is described in this section. 

No claim is made for optimality or near optimality in this approach. 
After many computations of the magnetic induction had been performed 
using this procedure, it became apparent that simplifications could be 
made, while still obtaining sufficient accuracy. However, these simp- 
lifications have not been tested, and the original approach to quadra- 
ture of the integral in (5) will be described. 

Measurements of the radial magnetic induction near the coil side 
face of the target in the EIDI prototype system provide some insight 
into a suitable numerical procedure for evaluating the integral in 
(1). The Fourier transforms of these observed fields are smooth 
(higher order derivatives with respect to frequency are small) with 
most of the energy confined to relatively low frequencies. This 
suggested the use of an IDFT to numerically perform the inverse Four- 
ier transformation in (1). Further impetus for use of the IDFT is 
given by noting that the DFT procedure for calculating the Fourier 
transform He O) of the current, described in the preceding section, 
yielded transform values at equally spaced frequencies. These are the 
appropriate frequencies for calculating a 1024 point IDFT, with a 5 
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microsecond sample time, of the Fourier space radial magnetic induc- 
tion. Hovever, numerical evaluation of the inverse Hankel transform 
in (5) at the 513 frequencies needed for calculating a 1024 point IDFT 
yielding a real sequence is prohibitively expensive. This expense was 
avoided by using the previously noted fact that the expected Fourier 
transform of the magnetic induction is smooth, vith most of the energy 
confined to lov frequencies. Accordingly, the inverse Hankel trans- 
form in (5) vas evaluated at only 77 frequencies, chosen to provide a 
reasonable representation of the behavior of the Fourier transform of 
the knovn radial magnetic induction. Spline function approximations 
vere then fitted to the real and imaginary parts of the calculated 
inverse Hankel transforms. The IDFT routine to calculate the desired 
time domain magnetic induction uses these spline functions to form 
B, (0, r, ) at the 1024 frequencies needed. 

With a procedure for evaluating the inverse Fourier integral in 
(1) available, a search for a suitable method for quadrature of the 
inverse Hankel integral vas initiated. The infinite upper limit of 
this integral vas simply replaced vith a suitably large finite upper 
limit, due to the controlling factor in the leading behavior of the 

— \q 

asymptotic expansion of the integrand being e. d . With the problem 
of the infinite upper limit gone, the remaining difficulties may be 
divided into tvo somevhat overlapping classes. 

Selection of a suitable algorithm for performing the quadrature 
is essential if accurate results are to be obtained. Most of the 
common quadrature algorithms are incapable of dealing vith oscillatory 
integrals vithout some help. Even vith the best of help, examples of 
problems vhere these algorithms completely fail abound. This has 
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given rise to highly specialized techniques for quadrature involving 
oscillatory integrands [27], [28], [29], (303. The use of one of 
these specialized algorithms vas avoided by writing a double precis- 
ion Gaussian quadrature routine [211, (311. This routine was then 
used to perform the inverse Hankel integration in (5) over the range 
0 < X < 10, then over the range 10 < X < 20, etc., stopping when the 
desired upper limit of integration had been reached. The results of 
these integrations, which can be interpreted as terms in a sequence 
that are to be summed, were then added together first using straight- 
forward addition, and then using the Euler series summation converg- 
ence acceleration algorithm DTEUL from the NASA Levis Research Analy- 
sis Center Software Library [343. For each radius at which the radial 
magnetic induction was evaluated, and for each of the 77 Fourier 
frequencies, good agreement was achieved between the results of these 
two summation methods. 

Several different finite upper limits were used to take the place 
of the infinite upper limit in the inverse Hankel quadrature. It was 
experimentally determined that upper limits greater than 1000 resulted 
in very little change in the calculated induction fields. 

The second problem area concerned the precision of the FORTRAN 
implemented on the 370. Single precision floating point word length 
on the IBM 370 is approximately 7 decimal digits (24 bit mantissa) 
(253. This is insufficient for nearly all complex scientific calcula- 
tions (323, (333. Double precision floating point word length on the 
370 is approximately 15 decimal digits (56 bit mantissa) [253, and vas 
used for all computations involving the inverse Hankel integral. 
Without the use of double precision arithmetic, inverse Hankel quadra- 
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ture was impossible due to smearing. 

Although the use of double precision greatly reduced the effects 
of smearing, it created problems of its own. A double precision 
Bessel function, and a double precision algorithm for K , ( K ) become 
necessary to evaluate (5>. Simply converting a single precision 
algorithm for K y ( X ) to double precision does not yield sufficient 
accuracy to obtain good results. Double precision Bessel function 
algorithms were unavailable at The Wichita State University. Double 
precision algorithms were written for J#(x) and J, (x). The double 
precision algorithm that was written for K / (X) is described in Appen- 
dix C. Finally, the complete program for calculating the time domain 
radial magnetic induction, given by (5), may be found in Appendix D. 
Figure 61 provides a comparison between the measured radial magnetic 
induction close to the coil side face of the target, and the predicted 
radial magnetic induction on the coil side face of the target (output 
from the program in Appendix D). 

6. 4. 3 Axial Magnetic Induction on the Coil Side Face of the Target 
The magnetic induction B, on the target next to the coil is 
calculated by performing an inverse Fourier transform on (22) in 
Section 4. 4 of Chapter 4. Then 

6 , (®, '{ 

r OO ~ 

^ 5 v (v-\ ixj- 
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where V(o, X , oj ) is given by (16) in Section 4.3.1 of Chapter 4. Sub- 
stituting (16) from Chapter 4 into the above. 
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The similarity between (5) and (6) is evident, even though the order 
of the inverse Hankel transforms is different. Because of this simi- 
larity, the algorithms developed to perform the quadrature in (5) were 
used, with only evident minor changes necessary, to perform the quad- 
rature in (6). A comparison at several different radii between the 
measured axial magnetic induction close to the coil side face of the 
target and the calculated axial magnetic induction on the surface of 
the target closest to the coil is shown in Figure 6J. 

6.4.4 Magnetic Induction on the Opposite Face of the Target 

Although it was not needed in the calculation of the force be- 
tween the coil and target, the magnetic induction on the side of the 
target away from the coil was calculated for comparison with the 
measured induction close to the same target surface. Section 4. 3 in 
Chapter 4 derived expressions for the total voltage V, and total 
current If at the junction between the metal line and air line oppo- 
site the current sources. These quantities correspond respectively to 
the azimuthal electric intensity and the radial magnetic induction on 
the surface of the target opposite the coil. 
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FIGURE 61 

B, VS. TIME, NEAR SIDE 
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Inverse Fourier-Hankel transformation similar to that in (5), but 
with the integrand corresponding to the current I, given in Section 
4. 3, yields the radial magnetic induction. Only minor changes are 
necessary in the program whose listing appears in Appendix D to calcu- 
late the induction. Figure 6K provides a comparison between the 
measured and calculated radial magnetic induction close to and on the 
target surface away from the coil. 

Inverse Fourier-Hankel transformation similar to that in (6), but 
with the integrand corresponding to the voltage V, given in Section 
4.3, yields the axial magnetic induction. Figure 6L shows the compar- 
ison between the measured and calculated axial induction close to and 
on the target surface away from the coil. 

6-5 FORCE VERSUS TIME 

It was shown in Section 4. 5 of Chapter 4 that the total force 
between the coil and target is given approximately by 

oo 

[&£(*, 0 ±) - r ( ^ 

o 

For use in this integral, radial magnetic induction was calculated at 
5 microsecond intervals at radii .01 inch, .2 inch, and in .2 inch 

increments up to a maximum of 2. 0 inch, on the surface of the target 

closest to the coil using the algorithms described in Section 6.4.2. 
Axial magnetic induction was also calculated at 5 microsecond inter- 
vals at radii .01 inch, .1 inch, and in .2 inch increments up to 1.9 

inch, on the same target surface using the algorithms described in 

Section 6.4.3. Cubic spline functions were then fitted to the squares 
of the radial variation of these two magnetic fields at times 0, 50, 
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100, 150, • . . , 900, and 1125, 1400, and 2000 microseconds for use in 
performing the quadrature indicated in (7), using an upper limit of 
1.9 inches. The quadrature was done exactly (to the limit of the 
machine arithmetic) on the spline function approximations to the 
squared fields by integrating the cubic polynomial form of the spline 
function between knots, and using the spline function coefficients in 
the result. Appendix E contains a listing of the FORTRAN program for 
performing this force calculation at t=50 microseconds. 

Figure 6M shows the force versus time calculated using this 
procedure. 


6.6 IMPULSE 


Equation (25) in Chapter 4, repeated below, was used to calculate 

P - ^ 


o 

the impulse delivered to the target by the coil. A spline function 
was fitted, with time as the variable, to the force calculated in 
Section 6. 5 above. An exact integration was performed on the cubic 
polynomials of the spline function between knots, with the result that 
the impulse calculated was .008 lb-sec. This is lower than the .012 
lb-sec impulse measured using a ballistic pendulum. It should be 
noted that the quadrature of the force integral in (7) above used an 
upper limit of r=1.9 inches instead of infinity. Although the induc- 
tion decays as the radius tends to infinity, an unknown part of the 
force has been ignored by not taking the induction fields at radii 
greater than 1.9 inches into account. 
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CHAPTER SEVEN 


CONCLUSION AND RECOMMENDATIONS FOR FURTHER WORK 

7. 1 CONCLUSION 

7.1.1 Summary 

A method of modeling the electrical system aspects of a coil and 

I 

metal target configuration resembling a practical Electro-Impulse De- 
Icing installation, and a simple circuit for providing energy to the 
coil, was presented. The model was developed in sufficient theoreti- 
cal detail to allow the generation of computer algorithms for the 
current in the coil, the magnetic induction on both surfaces of the 
target, the force between the coil and target, and the impulse deliv- 
ered to the target. These algorithms were applied to a specific 
prototype EIDI test system for which the current, magnetic fields near 
the target surfaces, and impulse had previously been measured. 

Coil impedance was the first quantity calculated using the algo- 
rithms. Agreement between the impedance calculated and the impedance 
measured was seen to be very good for the resistive part, and reason- 
able for the reactive part. Despite the simple model used for the 
circuit providing energy to the coil, excellent agreement was obtained 
between the predicted and measured coil current. 

Measured and predicted magnetic induction fields were not di- 
rectly compared, due to the fields having been measured close to, but 
not on, the target surfaces using spatial averaging methods. The only 
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fields calculated vere on the target surfaces, for use in calculation 
of the force between the coil and target. Nevertheless, there was 
seen to be very reasonable agreement between measured and calculated 
magnetic fields. The character of the time variation of these fields 
changed considerably with radial distance from the axis, and the 
algorithms correctly predicted these changes. 

Calculation of the impulse easily provided the greatest disagree 
ment between a predicted and measured quantity, with a -33X error in 
the calculated impulse. Impulse was measured using a ballistic pendu 
lum containing the metal target, so, for this measurement, the metal 
target was no longer held rigid when the capacitor was discharged 
through the coil. This does not satisfy the model assumption of a 
stationary system. However, the period of the pendulum was suffic- 
iently long that, during the time interval when most of the force was 
developed on the target, negligible motion of the pendulum should 
theoretically have occured. Motion of the target during the measure- 
ment of the impulse is not felt to be a satisfactory explanation for 
the discrepancy between the measured and calculated impulse. It was 
mentioned in Section 6.6, where the calculation of the impulse was 
described, that the infinite upper limit in the integral yielding the 
theoretical impulse had been replaced with a finite upper limit for 
the purpose of quadrature, thus incurring an unknovn error. Our 
procedure for quadrature of the impulse integral is felt to be the 
most likely source of error in the calculated impulse. 

7.1.2 Contributions by the Author 

Reference C41 provided most of the basic methods used in this 
dissertation in modeling the interaction between the coil and target. 
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It was shown that an error had occured in the definition given in C4] 
of the characteristic impedance of a Hankel space transmission line, 
and the corrected definition was used in the theoretical development 
of the model. For many transmission line calculations, this error 
does not become evident because impedances occur in ratios (e.g., the 
voltage and current reflection coefficient calculations in Section 
4.3.1). 

Reference 14] provided no method of calculating the axial magnet- 
ic induction B, <z, r, t) from the transmission line model. An integral 
solution for this field, in terms of the transmission line voltage, 
was derived in Section 4. 4. This solution allowed the use of nearly 
all of the algorithms developed for the calculation of the radial 
magnetic induction B, (z, r, t) in calculating the axial magnetic induc- 
tion B. <z, r, t). 

The most significant contribution of this work was the develop- 
ment of FORTRAN algorithms for performing the inverse Fourier-Hankel 
transformations yielding the induction fields, described in Section 
6.4. While the author made no theoretical contributions in developing 
these algorithms, several diverse results from the fields of numerical 
analysis and approximation theory had to be applied in concert to 
create a working algorithm. During this development, an algorithm for 
the generation of Struve functions of the first and second orders was 
developed that, according to a computer search of the literature, is 
the most accurate reported. 

7.2 RECOMMENDATIONS FOR FUTURE WORK 


If the methods developed in this dissertation are to be applic- 
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able as design tools for EIDI systems, the algorithms for the calcula- 
tion of the magnetic induction fields should be made more efficient. 
These algorithms, listed in Appendix D, take nearly half an hour of 
CPU time on the IBM 370 to calculate either the axial or the radial 
magnetic induction versus time at one specific radius. Both the 
radial and axial components are required, each at eleven different 
radii, for the force calculation described in Section 6.5. A minimum 
of eleven hours of CPU time is too great for the evaluation of the 
force-time profile of a proposed EIDI configuration. The least desir- 
able feature of the methods presented in this dissertation is the 
inordinate CPU time required for the calculation of the induction 
fields. A sophisticated convergence acceleration routine could per- 
haps be devised specifically for more economic calculation of the 
induction fields. 

Human intervention is required to proceed along the Analysis Flow 
Diagram of Figure 1 (page 4) in evaluating a specific EIDI system. 

This is because the outputs from the various programs, represented by 
such Figure 1 quantities as the coil impedance Z (co ) , the current i(t> 
and I(aO), are not in the form required as the inputs for the program 
at the next block in Figure 1. A significant amount of design automa- 
tion could be accomplished by vriting one long program, using as a 
skeleton the programs developed for this dissertation, that will take 
as input the geometry of the coil-target configuration and the circuit 
used to provide energy to the coil, and provide as output a force-time 
profile and the total impulse delivered to the target. 

Not all possibilities have been exhausted in the search for an 
analytical (or semianalytical ) solution to the EIDI design problem. 
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Perturbation techniques have been suggested as a possible method for 
analytical Hankel inversion in the calculation of the magnetic induc- 
tion , and should be investigated, as this is the area suffering the 
greatest computational expense. Furthermore, analytical solutions (if 
sufficiently simple) are often capable of providing insight into a 
problem not provided by numerical solutions. 



oo o oo o o r> o 


original rage is 

of poor Quality 


FORTRAN SOURCE CODE FOR CALCULATING 
COIL RESISTANCE INCREASE DUE TO TARGET 


C FILE: FFS21R FORTRAN 

C PR OCR AH TO ESTIMATE RCOIL, USING VALUES OF LAMBDA FROM 0 TO 680. 

C 

C VALUES OF THE E-FIELO ARE GENERATED DIRECTLY USING THE COOE 
C FROM THE FILE FEB2 FORTRAN. VALUES OF K.PRIME ARE OBTAINED FROM A 

C SPLINE FUNCTION APPROXIMATION. VALUES OF RCOIL CALCULATED WILL 8E 

C USFO IN A NUMERICAL INVERSE FOURIER TRANSFORM TO CALCULATE CIRCUIT 
C QUANTITIES. 

C 

C INITIALIZE 
C 

CCMMCN W 

DATA PI / 3 . 141 592 65/ 

TwCP 1*2 . *P I 

MAIN PROGRAM USING 0UANC8 


A=0 • 0 1 
U=680. 

PRINT OUTPUT HEADINGS 
WR[TE(6i50I A ,U 

50 FORMAT!' RESULTS OF NUMERICAL INTEGRATION FROM LAMBDA *',F4.1,« TO 
l'»F6.1/I 

GENERATE FREQUENCIES IN HERTZ 
DC 10 1*1,4 
OC 10 J=1 ,4 

GFNERATE FREQUENCIES INAl-2-4-7 SEQUENCE 
F=10.**I 

IF! J.EQ.2) F = 2. * F 
IF ( J . EC.3 ) F*4. *F 
TFfJ.EQ.4I F=7.*F 
W=TWCPI*F 
FFLERR=l .E-8 
A8S ERR = 0 . 

CALL CUANC8 I A, U, A BS ERR, RE LERR ,RE SULT , ERREST , NOFUN , FI AG ) 

INTEGRATION OONE - OUTPUT COIL RESISTANCE CALCULATED 
RCOI L=T WOP I* RESULT* 1000. 

10 WPITE!6, 1010 I F.RCOI L, ERREST 

101C FORMAT!' FREQUENCY* * ,F6.0,5X,' RCOIL* • , F 7. 3, • M ILL IOHMS ' , 5X, 
I'EPPEST* ' , E 12. 5 ) 

STOP 

END 

COPIED SL3R0UT I NE QUANC8 FOLLOWS 

S'JBRCL'T INE QUANC8IA, B , AB SERR , REL ERR , RESULT , ER REST , NOFUN, FL AG) 

REAL FUN, A, 8, ABSERR, RELERR, RESULT, ERREST, FLAG 
INTEGER NOFUN 

REAL V,0,Wl ,W2,W3 ,W4, AP F A , XO , F 0, S T ONE , ST E P , COR l 1 , TEMP 

REAL CPREV,QN0W,Q0IFF,QLEFT,ESTERR,10LERR 

REAL CRIGHT ! 3 1 »,F(I6),XII6),FSAVE(8,30),XSAVEI8,3P) 

I NT F G c R levmin,ievmax»levqut,ncmax, nuf in.lev.nim, i, j 

lev-in = 1 

LEVMAX = 30 

LEVCUT = 6 

NC“AX * 5000 
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NQF IN = NCMAX - 8 * I L E VM AX-L E VOUT ♦ 2* * I LE VO UT* 1 J ) 

WO = 3956.0/14175.0 
Wl = 23552.0/14175.0 
W2 = -3712. 0/14175.0 
W3 = 41984.0/14175.0 
. W 4 = -18160.0/14175.0 
FLAG = 0.0 
P.FSULT = 0.0 
CTPll = 0.0 
ERREST = 0.0 
AREA = 0.0 
NC FUN = 0 
1FIA.FC.B) RETURN 
L F V = 0 
N I u = l 
XO = A 
X (16) = 8 
OPREV = 0.0 
FO = FUN(XO) 

STONE ■ (B-AI/16.0 
X ( 8 ) = (XO ♦ X( 16) 1/2.0 
X ( 4 ) = (XO ♦ X ( 8 ) 1/2.0 
X ( 1 2 ) = ( X ( 8 ) + X ( 16 ) J/2.0 

X ( 2 ) * (XO ♦ X ( 4 > ) / 2 . 0 
X(fc)-(X(4)*X(8))/2.0 
X(10) = ( XI8) ♦ X( 12) 1/2.0 

XI 14) = ( X I 1 2 ) ♦ XI1HI/2.0 
OC 25 J = 2,16,2 

c l J) = FUN ( X ( J ) ) 

25 CONTINUE 
NOF'JN = 9 

30 X ( 1 ) = 1X0 ♦ X(2) J/2.0 
F ( l ) = F UNI X ( 1 ) ) 

CO 35 J = 3,15,2 

XCJ) = (XCJ-Il + X(J*l))/2.0 
FIJI = FUNIX ( J ) ) 

35 CENT 1 NUE 

NOFUN = NOFUN ♦ 8 

STEP = I X ( 1 6 ) - X0I/16.0 

3LEFT = ( WO* ( F 0 + F(8)l ♦ W 1 * ( F ( 1 ) *8(7)1 ♦ W2*(F(2) ♦ F ( 6 J » 
l ♦ W3* ( F l 3 ) ♦ F(5)) + W4*F(4) ) *$TEP 
QP IGHTI LE V* 1 ) = ( WO* (F t 8) *F( 16 ) )*W1*( F(9)*F ( 15 ) ) + W2*(F ( 10) *F( 14) ) 
1 ♦W3*(F(ll)+F(13))*W4*F(l2))*STEP 
CNOW = CLEFT ♦ 0PIGHT(LEV*1) 

CO IFF a QNOVi - QPREV 
AP = A = AREA ♦ CDIFF 
ESTERR = ABSIQC IFF J/1023.0 

TCLERP = AM AX 1 ( A3SFPR ,RELERR*ABS I AREA) )*( STEP/STONF ) 

IFIlEV.LT .LEVNIN) go TO 50 
IFCLEv.GE .LCVM/.X I GO To 62 
IFINOFUN.GT.NOF IN) GC TC 60 
IFIESTERP.LF.TOLERR) GO TO 70 
50 NI** = 

LEV = LEV +1 
00421=1,8 

FSAVEt I, LEV) = F( 1+8) 

XSAVE ( I, LEVI = X( 1+8 ) 

52 CGNT INUE 

OPREV * OLEFT 
00 55 I = 1, 8 
J = -I 

F ( 2*J* 18 ) = F ( J+ 9) 

X (2* J*18 ) = X(J+9) 

55 CONTINUE 
GO TO 30 

60 NOF 1 N = 2*N0F IN 
LEV*IAX = LEVOUT 
FLAG=FLAG*( B-XO) /(8-A ) 

GO TC 70 
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62 FLAG=FLAS*1.0 
70 RESULT=P£SULT*0N0w 
£RRFST = EKREST«-ESTtRR 
CCR11=COR11+QDIFF/1023.0 
72 1FININ.EQ.2MNIM/2) ) GO TO 75 
NIM= MM/ 2 
LFV = LEV-l 
GO TC 72 

75 N I M = MM ♦ l 

IF(LEV.LE.O) GO TO 80 
CPREV = QR ! GHT (LEV) 

xo = x( it ) 

FO = F ( 16 > 

DO 78 I = 1, 8 

F ( 2* I ) = FSAVEI l.LEV) 

X(?*I) = XSAVE ( l.LEV) 

78 CONTINUE 
GO TC 30 

80 RESULT = RESULT ♦ CORll 
IFIERREST . EQ.O. 0 ) RETURN 
82 TEMP = ABS(PESULT) ♦ ERREST 

IFITEMP.fJE.ABSIf ESULT ) J RETURN 
ERREST = 2.0*EPPEST 
GC TC 82 
END 
C 

C COPIFD SUBPROGRAM SEVAL FOLLOWS 

REAL FUNCTION S F VAL I N , U . X , Y , B ,C t D ) 

REAL L,X(N),Y(N) ,B(N) tC(N),0(N) 

DATA 1/1/ 

IF(l.GE.N) 1=1 
JFIU.lT.xmi GOTO 10 
IFIU.LE.XII + U ) GOTO 30 
10 1 = 1 
J = N+1 

20 K = t I ♦ J ) / 2 

IF (U.IT.X(K) 1 J - K 
IFIU.CE.X ( K J ) I=K 
I FI J.GT. I *1 ) GOTO 20 
30 DX=U-X< I I 

SEVAL* YC n + DX*IB( I)*OX*( C( I) *OX*Dt m ( 

RETURN 
END 
C 

C SUBPROGRAM FUN(LAMB) TG EVALUATE THE INTEGRANO FOR QUANC8 
FUNCTION FUNILAM8) 

REAL LAMBtKPRI ,MUC/12 .5 663 7E-7/ , L AMBDAI 48>,KPRIHE(48),B(4*),C(4ei, 
10(48) 

COMMON W 

COMPLEX CMPL X ?C SORT ,CEX P , Z PDZ A, I N , ZXDZA , R HOT , H, E t T1 , HP 2 , HPO.CC t CG t 
1 BC t BG t AC i AGt ZA.ZASQ, E2t El t E3 t H3, E PO, HR1 
DATA SIGMA/1. 74E7/,H03?. 00159/, G/ .00278/ , 00/ . 0008 128/ 1 IFLAG/l/ 

C REAC IN SPLINE INTERPOLATION DATA FROM FILE KPSPL C QEFF 
C FIRST CHECK TO SEE IF COEFFICIENTS ALREADY READ 
I F ( IFLAG.NE.l ) GOTO 10 
00 1 1=1,48 

l READ t 2, 10C0 ) LAMBOAI I ) .KPRIME tl) ,B(I) ,C( I ) ,0( I) 

1000 F0RMATIF6.1,4(1X,E15.6I ) 

C FIRST, EVALUATE THE PEAL ANO IMAGINARY PARTS OF THE E-FIElDS REOUIRFO 
C FOP THE INTEGRATION - USES CODE FROM FILE FEB2 FORTRAN 
C FOLLOW STEPS OUTLINEO IN "METHOD" 

IF(LA«B.EO.O.) LAMB* l . E-20 

10 ZPDZA=( l. ,0. I/CSOPT ( CMPLXI l. , W*MUO* SI GMA/ ( L AM B*L AMB ) ) ) 
T=EXP(-LAMB*H03) 

I N=C MPl X ( I T*T*T*1 ) *EXP I -LAMB *G> ,0. ) 
RHOT=((I.,0.)-ZPDZA)/((1.,0.)+ZPDZA) 
Tl=RHCT*LEXP(CMPLX(-2.*LAMB*DD,0. I/ZPDZA) 

zxczA=m .,o. i+ti )*zpoza/((i . ,0. i-ti j 

H3 = I N / 1 ( 1 .,0. )*ZXOZAI 
ZA=CM°LX(0., W*HUO /L A MB I 
E3=H3*ZX0ZA»ZA 
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C CALCULATE EP2 

T*FXP(LAMB»GJ 

AG=C MPi_X( ( T* l./T) n. , 0. I 

9G=Cf PLXI IT-l./T »/2. ,0. )*2A 

zaso=za*z a 
CG*BG/ZASC 
E 2 = A G*E 3*BG*H3 
HR2 = CC*E3-*AG*H3 
EP?R=BEAL ( E2 I 
EP2 I* A IMAG (£2 ) 

C CALCULATE EP1 

T=EXP(LAMB*HD3) 

AC=CPPLX C (T+1./T1/2. ,0. ) 

8C=CMPLX< (T-l./T ) /2. >0. )«ZA 
CC*BC/ZASG 

f l=AC*£2»aC*(HP2- ( 1. .0. n 
HR1*CC*E2*AC* ( HP 2- ( l ..0.1) 

EPin=PEALtEl) 

EP1 I* A IMAG (El I 
C CALCULATE EPO 

EPO* AC* El *BC*( HP 1- ( l . ,0. ) 1 
EPOfi=FEAL (EPO) 

EP01=AIMAG(EP0) 

C EVALUATE KPRI USING T HF SPLINE INTERPOLATION 
KPR!=SEVAL (48, LAMC, LAMBDA, KPRIME , B,C,D) 

C NOW EVALUATE THE INTEC.RANO AND PFTURN 
FUN=LAMB*KPP I *KPR1 *( E P2P ♦EP1 P ♦EPOB ) 

I p LAG=2 

RETURN 

END 


RESULTS OF NUMERICAL INTEGRATION FPOM LAMBDA = 0.0 TO 680.0 

EPP EST= 0.960R5F-13 
ERREST= 0. 1 1895E-1 1 
EPPE ST* 0.14B12F-11 
ERR E ST* 0. 1 6755E-10 
ERREST= 0.23993E-10 
ERRE ST= 0. 30565E-09 
ERR E ST = 0.3R682F-09 

ERREST = 0.34003E-09 

ERREST* 0. 33792E- 09 
ERREST * 0.52783F-08 

ERREST* 0. 91 543 E — C 8 
ERREST* 0 . 58640E-08 
ERREST* 0.64146E-OP 
EPPEST* 0.39609E-08 
ERREST* 0.49584E-08 
EPPEST* 0.64945E-08 


FREQUENCY* 10. RCOIL= 0.004 MI LLI OHMS 
FREQUENCY* 20. PCOIL= 0.017 MILLIOHMS 
FREQUENCY* 40. RCOIL* 0.068 MI LL I OHMS 
FREQUENCY* 70. RCOIL* 0.207 MI LL I OHMS 
FREQUENCY* 100. RCOIL* 0.418 M I LL I OHMS 
FREQUENCY* 200. RCOIL= 1 . 589 MI LL I OHMS 
FREQUENCY* 400. RCOIL* 5.592 MILLIOHMS 
FREQUENCY* 700. RCOIL* 13.341 MI LLI OHMS 
FREQUENCY* 1000. RCOIL* 22.783 MILLIOHMS 
FREQUENCY* 2000. RCOIL* 47. 3 73 MI LL I OHMS 
FREQUENCY* 4000. RCOIL* 70.588 MILLIOHMS 
FREQUENCY* 7000. PCOIL* 81.442 MILLIOHMS 
FREQUENCY* 10000. RCOIL* 85.525 MILLIOHMS 
FREQUENCY* 20000. PCOIL* 92.622 MILLIOHMS 
FREQUENCY* 40000. RCOIL* 109.270 MILLIOHMS 
FPEQUENCY* 70000. PCOIL* 142.962 MILLIOHMS 
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APPhXDIX ij 


C 


FORTRAN SOURCE CODE FOR CALCULATING 
L CURRENT BEFORE CLAMP DIODE CONDUCTION 


PROGRAM TO PERFORM NUMERICAL FOURIER INTEGRAL INVERSION TO CALCULATE 
C THE CURRENT IN THE COIL IN THE WING CIRCUIT (BEFORE THE 

C DIODE CLAMP CUTS IN). THIS PROGRAM INCLUDES THE FREQUENCY 

C DEPENDENCE OF THE COIL RESISTANCE AND INDUCTANCE. OUTPUT FROM THIS 
C PROGRAM WILL BF USED TO CALCULATE THE SPECTRUM OF THE CURRENT. 

C THIS PROGRAM IS ADAPTED FROM SERCUR FORTRAN. 

C 

C 0 MO L EX F1(2048),CM?LX,C0NJ3 

R C AL L.LSIMID r , 

DATA N/2D48/, T/5.E-6/,C/fcOO.E-6/,RCA3LE/. 05*/ ,RDC /. 047 /, L SI ft I D/23 . 
*E-6/, PI/3. 141593/ 

TTOT=FLOAT(N)*T 
D EL TAF = I . /TTOT 
DELTAW=2.*PI*DELTAF 
K ODE = -1 

FORM COMPLEX SPFCTFUV MATRIX 


I END=N/2+l 
Dn 5 1 = 1 , IEND 
W=FLOAT( I-1)*DELTAW 
C CHECK FOR FREQUENCY w=D 
IFtw.NE.O.) GOTO 3 
R=RCABLE+RDC+.004E-3 
L=LS IMID+17. 865E-6 
GOTO A 

C W .NF. 0 - CALCULATE RESISTANCE R AS THE SUM OF 3 TERMS: A CONSTANT 

C D.C. TERM (RCABLE); SKIN EFFECT IN COIL (RAC); AND LOSS IN SKIN OF 
C. WING (COUP) 

3 R=PCABLF+RAC (W)+C3ILP (W) 

C CAlCULATf INDUCTANCE AS THE SUM OF 2 TERMS: A CONSTANT TERM DUF TO 
C the S I m.mqnds AND IDLER (LSIMID) ; AND A TERM DUE TO THE COIL IN THF 
C WING (Cf'ILl) 

L=LSIMID+COILL (w) 

C NOW FORM THE CONSTANTS USED IN EVALUAT ION OF THE FOURIER TRANSFORM 
C OF THE TIME VARYING PART OF THE CIRCUIT CURRENT 

A P1=L*C 
P2=R*C 

C NOW FORM THE MATRIX VALUES 

5 F 1( 1 )=( 1. ,0. ) /CMPLXl l.-W*W*Pl »w*P2) 

I STAR T=f,/ 2+2 
DO 7 1= I S TAP T , N 
7 FI ( I )=CUNJG( FI (N+2-* ) ) 

C 

CALL FFT(KOOE,N, DELTAF.FI, 13) 

C EXPRESS TTCT IN MRLISFC, T IN MICROSEC 
TT0TP=TT0T*'1C0D. 

TP= t *I.E6 

C BEGIN OUTPUT OF RESULTS 

W D ITF(t»20I) N,TTOTP,TP 

?01 FlrvaT (• 1 NUMBER OF P 3 INTS = • , I A/ • TOTAL TIME I NTERVAL= *,F7.?,* MR 
1 L I SEC ’ » / ’ SAMPLING PERIOD= ',F7.3,' MICROSEC'//' OUTPUT FROM OROGR 


IDD 


20 
2 on 
10 


! A m M^Y22A FORTP AN* II ) 

WPI TE ( 6 , 100) 

FORMAT!* CIRCUIT CURRENT 
1 IME(MStC) ' .5X , • TOTAL CURPFNT 
FSC=400.*C 
DO 20 1=1 ,100 
DT=FLCAT( I-l.)*T*1000. 

FI ( I )=CMPLX(FSC,0.)*F1( I ) 
WR!TE(6,200) DT,F1(I) 

FORMAT (IX ,F12.3,1X,2(F15.5,5X)) 

STOP 

END 


FREQUENCY DEPENDENT L, 
RESPONSE* I ,T23 ,'REAL* , 15X 


P'//* T 
* 1MAG’/) 
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FUNCTION R AC ( K ) 

SUBPROGRAM TO EVALUATE SKIN EFFECT RESISTANCE RAC OF COIL AS A 
FUNCTION OF THE RADIAN FREQUENCY W 

COM.RLFX CMPL X.CEXP , TAU, TAUHD2 ,T1 ,T 1R 
PEAL MU/1.256637E-6/ 

DATA H/.00477/,SIGM4/3.48E7/»R0C/.047/ 

OELT A= SORT (2 . / ( W*MU^S IGMA } ) 

TAU=( I., 1.) /CMPL Xf DELTA, 0.) 

T AUHD2=TAU*CMPLX ( H/2 . ,0. ) 

T1=CEXP( TAUHD2) 

Tl-=(l.*0.)/Tl 

R AC- R DL*P EAL ( TAUHD2M TI +T 1 P > / ( T 1 -T IP I) 

RETURN 
e i\p 

FUNCTION COILR(W) 


SUBPROGRAM TO EVALUATE COIL RESISTANCE OOILR (DJE TO PROXIMITY OF 
WING SKIN) AS A FUNCTION OF THE RADIAN FREQUENCY W 


PEAL FREQ(24),pr.0rL(24),BR(24),CR(24) ,DR(24) 

DATA P1/3.1415 C 3/,IFLAG/1/ 

IF (IFLAG.NE.I) GOTO 5 

READ IN COEFFICIENTS FROM FILE P CO I L COEFF FOP SPLINE FUNCTION 
INTERPOLATION 
DO 1 1=1,24 

1 READ! 2» 100) FRE 0 (1) , RCO I L II) , OR (I) ,C R (I) , DR ( I ) 

105 FORMAT ( F 8 • 0 , A ( IX, E15.6) ) 

T WO I = 2 . * P I 
IFLAG=-1 
5 F = 'w / T wO P I 

IF (F.GT.10.) GOTO ! C 
C Cl LR = • 004E-3 
RFTURN 

10 COIL R=SEVAL( 24 , F ,FPfQ , PCOIL.BR, CR ,DR I. E-3 
RETURN 


F NIP 

PEAL FUNCTION SE V AL l N ,U . X , Y , B ,C , D > 


COPIED SUBPROGRAM SEVAl FOLLOWS - EVALUATES SPLINE 

P R A l U,X(N),Y(N) .BIN) ,C(N) ,D< N) 

data 1/1/ 

IF(I.GE.N) 1 = 1 
IFIU.LT.XI I ) ) GOTO 10 
IFlU.LE.XII+1) ) GOTO 30 
10 1=1 
J = N+1 

20 K=(I+J)/2 

IFIU.LT.X (K) ) J=K 
IF(U.GE.XIK) ) 1 =K 

IFIJ.GT.I+l) GOTO 20 
30 OX=U-X(I) 

SEVAL =Y( I )+DX*( B( I ) +DX*( C ( I )+DX*DC I) ) ) 

RETURN 

END 

FUNCTION COILL(W) 


SUBPROGRAM TO EVALUATE COIL INDUCTANCE COILL AS A FUNCTION OF THE 
RADIAN FREQUENCY W 


PEAL FREQ (16) , LCOIL ( 1 6 ) , BL ( 16 1 , CL ( 16 ) ,DL( 16) 

DATA PI/3. 141593/, IFLAG/1/ 

IF (IFLAG.NE.I) GOTO 5 
C READ IN COEFFICIENTS FROM FILE LCOIL COEFF FOR SPLINE FUNCTION 
C INTERPOLATION 
DO 1 1=1, 16 

l REAO( 3, 100) FREQ(I), LC01LII) ,BL( 1 ) ,Cl( I ), DL( I ) 

IPG F0R v AT(r6.0,4< IX.E15.6)) 

TwOP I = 2 • v P I 
IFIAG=-1 
5 F = W/ T WOP I 

IF (F.3T. 10. ) GOTO 10 
C0ILL=17. 865E-6 
RETURN 

10 ir(F.GE.70.E3) GOTO 20 

CQILL=SEVAL(16,F,FREQ ,LCG I L , BL, CL , DL ) *1 .E-6 
RETURN 

20 COI LL=9.263E-6 
RETURN 
END 
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SUEPOLT I N E F F T ( KO DE , N , DE LT A , X , * ) 

PRCOPAV TO IMPLEMENT ThF FAST FOURIER TRANSFORM WHEN THF NUMRFR r>F 
DATA PTINTS tS AN 1 NTEGF R POWER OF TWO. TAKEN FROM PAGES 264-265 
Of "M r THOOS CF DISCRETE SIGNAL AND SYSTEM ANALYSIS", EY JONG. 


C 0 M P L FX X CM > ,W,X 1 ,CMPLX 
JP = 0 
N 1 = N 

5 N2=Nl/2 

IF(N?*2.NE.N1 ) GC TO 100 
t L ‘ = I R 
N l=N2 

IF (N1 .GT . 1 ) GOT 0 5 
PN = t. 2831 85/N 
L = N / 2 


I R 1= I P- 1 
K 1 = 0 

DO 30 I S= i , I Ft 
15 DO 20 1=1, L 
K=K 1+1 
KPL=K*L 

AM=KBITR(K1/2**IR1,IR) 
IFCAM.NE.O.I GOTO 18 
X 1 =X ( KPL ) 

GOTO 19 
13 AP.G = AM*PN 
C =C OS ( ARG ) 

S=-KODE*S IN ( ARG ) 
W=CMPLX(C,S) 

X 1 = W*X (KPL ) 

19 X ( KPL )= X { K ) -X 1 
X (K ) =X( K ) +X1 

20 K1=K1+1 
K l = Kl+L 

IF(Kl.LT.N) GOTO 15 
K 1 = 0 

I Rl= I Rl-1 
30 L=L/2 

DO 40 K= 1 , N 
K 1 =K3 I TR ( K— 1 , IRH-1 
IFIKl .LE.K) GOTO 40 
X 1 = X ( K ) 

X ( K ) = X ( K 1 1 
X IK1 J =X1 
40 CONTINUE 

IFIOELTA.EQ.l .) RETURN 
no 50 K=l,N 
50 X (K)=DELTA*X (K) 

P F TURN 

100 W R I T E ( 6, 101 ) N 

101 FORMAT!//, » *** N= ',16,' 
1 ' ) 

RFTURN 

END 

FUNCTION KBI TP ( K , I R ) 


K 3 1 TR = 0 

do” l 1=1, IP 
K2=K 1/2 

KB1TR=2*KBITR*K i-2*K2 
1 K 1 =K2 

. return 

F ND 


IS NOT A POWER OF 


FFT PUN ABORTED *** 



APPENDIX C 


CALCULATION OF K'< X ) 


The function we have called K 7 (X > arises from the transformation 


of real space current discs into Hankel space. By definition, 

20 


30 . . 

R 1 


(.) 


Straightforward quadrature of (1) using the Newton-Cotes algorithm 
QUANC8 was initially performed for generating values of K*< X) for use 
in approximating K'( X) with a cubic spline function. This provided 
sufficient accuracy for use in numerical calculation of coil imped- 
ance, as described in Section 6. 2. Attempts at calculating the radial 
magnetic induction using this spline function approximation for K , (X) 
were a total failure. This appendix describes a procedure for calcu- 
lating the function K' l X) that is accurate to 14 digits when imple- 
mented in double precision FORTRAN on an IBM 370. 

The integral in (1) can be evaluated in closed form in termB of 

named special functions. Change variables, 

Xx - ^ 4 dx = ^ /X 


* *, * 
- * 


% " 

^ «■ X 


so that (1) becomes 


\ O 

< M -- a~Z 


«■> 


f‘ * ^ * 

Xfc. 
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kVx'i = ^5 j aVa 1 h < s > 

A 1 "'xfc, 

From reference C18], 

where H 0 (x) and Hi (x) are Struve functions of orders 0 and 1 respec- 
tively C18], CIS], and J 0 (x) and J t (x) are Bessel functions of the 
first kind of orders 0 and 1 respectively. Using this result in (2) 
yields 

kW - 1>. W - W‘3'3 

In order to use this result, double precision algorithms for generat- 
ing the Bessel and Struve functions must be available. Double preci- 
sion Bessel functions are readily available, but Struve functions are 
not. A computer search of the literature resulted in a reference to a 
Naval Research Laboratory report that contained FORTRAN source code 
for generating integer order Struve functions with positive arguments 
[36]. We were unable to make use of this code because it used subrou- 
tines to which we did not have access. 

For several years, mathematicians have been avare of the desira- 
bility of using truncated Jacobi series of Chebyshev polynomials for 
numerical approximation of various special functions C37], C383. Luke 
C 37 3 provides coefficients b* and c„ for Chebyshev series expansions 

* t , 1*1 a 8 

• £^(f) , 


'"■A 

(*) 

3 ‘ U| 


l x l - & 
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where T a „<x) is a Chebyshev polynomial of the first kind of order 2n. 
All coefficients having magnitudes greater than 10'*° are given, 
allowing generation of H 0 <x> and H t (x) with 20 decimal digit accuracy 
for arguments whose magnitude is less than 8 [39]. Using the identity 
T a . (x) = T* ( 2x* -1 ) 

and retaining sufficient terms for 15 digit accuracy, these series 
expansions become 

4 2 , ix.i ; « (*) 

n-o u J ) 

Then for arguments X such that XR, <8, the expression on the right 
hand side of (3) is evaluated using (4) and (5). Chebyshev polyno- 
mials are evaluated in double precision arithmetic using the subrou- 
tine DCNPS from the NASA Lewis Research Analysis Center Software 
Library [403. 

Luke also lists coefficients d„ and e. for the series 


14, ' 

- \(A 

- SL 
~ TTX 

L ^ li 

:) , * 1 * 

U,M - 

Y i'x.'i 

. 3. 

~ t r 

L (l 

:) , 


for 15 digit accuracy. The functions Y»(x) and Y t (x) are Bessel func- 
tions of the second kind of orders 0 and 1 respectively. Rather than 
use these series directly to evaluate H 0 (x) and H t (x) for arguments 


greater than 8, some simplification is possible. Writing 




nr X. *-o *» 
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- nr 


(|) + Y,M = A + 


and substituting these expressions for H 0 (x) and H, (x) into the ex- 
pression J, (x)H 0 (x)-H, (x) Jo (x) appearing in (3), and suppressing the 
argument x, ve have 

J, Ho - H.Jo = J t C Ao + Yo 1 - CA, «-Y, ]J» 

s J| Ao — A| Jo * Jt Yo — Yt Jo 

= Jt Ao ~ At Jo * ' ~ 
nr x. 

using a veil known property of the Wronskian of Bessel functions of 
the first and second kinds C18]. Thus, the expression on the right 
hand side of (3) can be written as 

when it is to be evaluated at an argument y = XR. > 8. 

Using these results, K <X) can be evaluated from (3) for non- 
negative values of the argument to at least 14 digit accuracy. 
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ORiGi'NAL ‘IS 

OF POOR QUALITY APPENDIX D 


FORTH AN SOURCE CODE FOR CALCULATING 
RADIAL TIME DOMAIN MAGNETIC INDXT1CN 
CN THE COIL SIDE SURFACE OF THE TARGET 


riLf.: waue .2 fortran vi.s 

i ru calculate the -»aoial co mponent of rwi; magnetic induction 
at r m f. coil side pact of the skin, at a radial distance a inch 
F lOH THE CUTL AXIS. ADAPTED FROM SOURCE. AUG2 VI. 0 DROUGHT TO HSU 
FROM NASA LEWIS AUGUST AS. 

SSZS# CAUTION 

DARRAY SHOULD HAVE G I “6 NS I ON NARP.AY, AND DHINTL SHOULD HAVE 
DIMENSION NINTL. IN THE FUNCTIONS UFCT(I) AND DFUN(DX), OARRAY 
SHOULD HAVE DIMENSION NARAAY. NARAAY*1=NINTL 
rN AOOITION, the array ICKFLG SHOULD 8 E dimensioned as 
I E S F L '. ( 2»77,NARHAY>. 

OOSSS CAUTION 

VI. 1 - ADDED the ABILITY TO STORE AND PRINT THE FLAGS ASSOCIATED 

WITH EACH INDIVIDUAL GAUSSIAN DU A DRAT UR 6 . 

VI. 2 - CHANGED AH SOLUTE AND RELATIVE ERROR REQUESTS TO l.D-12 IN 

INVERSE HANKEL TRANSFORM QUADRATURE. 

CHANGED PRINT FORMATS TO ALLOW SPOOLING OUTPUT TO MY READER. 
VI. I - ADDED OUTER LOOPING A 3 1 L I T Y TO CALCULATE MAGNETIC INTENSITY 
AT SEVERAL RADII. 

IMPLICIT REAL-6 < 0 > 

COMPLEX CURA (77) ,CUBNT,CFMR( 10 Z 4 ) . C MUO/ ( 1 2 • 5& 6 37 E- T , 0 . > / 

REAL'S DARR AY( 100) , OHINTLI 10 l ) . DHRROI 77 ) , DH IROI 77) 

RE AL-4 DELTA 

real MONTH/’JAN «/«DAY/» 14 •/ t VERS/»V1.3*/ 

REAL OME G A ( 77 ) » HR R 0( 7 7 ) « HI RO ( 7 T ) » HR S B( 7 7 ) » HR SC ( 77 ) » 
'5MRSD(7T),hISH(T7),HISC(77),HISD(77) ,R INCH/. 01/ 

INTE'.FR I ERFLGI 2* 77 , 100 ) 

OaTA NARRAY/100/.NINTL/101/ 

EXTERNAL 0 F UN .DECT 

COMMON DARR A Y, CUPNT , IFUN»R » 7W0PI F 

BLOCK 1 - READ COIL CURRENT SPECTRUM (STEP T) (POSITIVE FREQUENCIES 

ONLY) FROM FILE CURKSPEC DATA (LOGICAL DEVICE 5). 
FREQUENCIES ARE STORED IN OMEGA { I ) WITH THE CORRESPONDING 
CURRENT SPECTRAL VALUE IN CURR(I). 

READ CURRENT SPECTRUM (HEADINGS HAVE BEEN REMOVED FROM THE FILE 
CUBA SPEC DATA) 
no to 1=1,77 

10 READ (S.IOOO) OMFGA( I ) ,CURR( I > 

1000 FQAMlT(TlR f El2.S»T3b,E13.S,TSl,El3.5) 

OUTER LOOP TO CONTROL RADIUS AT WHTCH MAGNETIC INDUCTION IS 
CALCULATED 


DO V» I ' AQUS=1S,20,2 

RINCH=FLOOT( IRADUS)/10. 

EXPRESS R IN METERS 
•>=R INCH-R.02S4 

"L 1C< 2 - CALCULATE INVERSE HANK EL TRANSFORM OF THE MAGNETIC 

INTENSITY. THIS IS DONfc AT EACH OF THF DIFFERENT 
FREQUENCIES QMEGA(I), 1 = 1.2.. ...77. THE FREQUENCY 
OMEGA IS PASSED TO THF SUHPROGRAM FUN THRU THE COMMON 
VARIABLE TWOPIF. the COMPLEX value of THE CURRENT 
SPECTRUM AT THIS FREQUENCY IS PASSED TO FUN THRU THE 
COMMON variarle CURNT. 

CALCULATE THE REAL** ARRAY OF UPPER AND LOWFR INTEGRATION LIMITS 

used in forming the sequence of magnetic intensities. 

0 S T E P = 1 0 • 0 0 
DO 1 1=1 .NINTL 

l DhINTL( I )=DFLOAT( I- l)*DSTF» 
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c print k-:j3in':,s 

WR I T i- ( *> * 40 )0 ) VE B S, "ONTH, DAY , DMl NTL( N I NTL ) 

R000 FiJRMAT(* 10UTPUT FROM WAUG..’ FORTRAN • I AU , 10 X , 84 • A A , / , • UPPER t.NVE 
SAjF HANK EL TRANSFORM I N T EC. A A T I Ofl LIMIT = *,013.3//) 

UR I T i. (*., R0J1 > R INCH, A INCH 

9001 fll’NIire FREQUENCY ARAL PART OF MAGNETIC INTENSITY* , T S 7 • ' I MAI', PA 
«ai of MAGNETIC INTF.NS I TV* , / , • ( HE R T 7 ) • « T 2 0 , • H ( A s • , F 3 . I , * , F ) * , T 3fl 

*, * I ER N*,Ts3,*H(R = *,FJ.l,«,F)*,T.<)l,*IFR N* > 

HE G I NNI NG OF MAIN LOOP TO CALCULATE INVERSE HANK EL TRANSFORM AT 77 
DIFFERENT FREQUENCIES 
DO 000 IMAIN=1,77 
THOPIF=OMEGA(IMAIN) 

CURNT=CURR( IMA IN) 

CALCULATE REAL PART OF HR(R,TWOPIF) 

I F UN = 1 

FORM THE SEQUFNCE OF PARTIAL INVERSE HANKEL TRANSFORMS BY INTEGRATING 
HETUEEN LAMBDA SU3 ( I ) AND LAMBDA SU3(I*1)» FOR 1=1 TO NARRAY (THESE 
LIMITS ARE IN THE ARRAY DM I N TL ) • 

SET tjp ABSOLUTE AND RELATIVE QUADRATURE ERROR REQUESTS 
D A [>, S f R - 1 . D - l 2 
DRtR-l*’)-l? 

INITIALISE TH ■ VARIABLE USED TO ACCUMULATE the results of the 
NUMERICAL INTEGRATION OVER EACH LAMBDA SUB-RANGE* 

DINTR=0.n0 
DO 2 I=l*NARRAY 

FORM THE LOWER AND UPPER INTEGRATION LIMITS 

oa=dhintl( I ) 

DB=DHINTL( 1*1) 

INTEGRATF ( T H r s E VALUES ARE PASSED TO THE FUNCTION DF C T THRU COMMON). 
RESULT OF TH = INTEGRATION IS STORED IN DARRAV(I) FOR FUTURE USE IN 
ThE EULE j CONVERGENCE ACCELERATION ALGORITHM DTEUL. 

CALL QGAuSS( DA ,DH, DFUN* DR£R* DArtSER* DARR AY( I)»DERR»IER) 

STORE ERROR FLAG 

IERFLG( I. IMAIN.I ) = I E R 

NOW SUM TH- INTEGRALS FOR C 0m® AR I SI ON WITH THE accelerated SUM TO 
BE CALCULATE! IN SUBROUTINE DTEUL. 

2 DINTS. rntNTR*DARRAY( t ) 

NOW USE THE f ul t R CONVERGENCE ACCELERATION ROUTINE TO BEST ESTIMATE 

THE REAL PART OF THE TOTAL INVERSE HANKEL TRANSFORM AT THE 
F'EJUCNCV TWOPIF. 

call oti;ul( dfc t, dsumr ,n4Rb ay , i . e - is , ierr »nr ) 

DHRRO( IMA IN) =OSUMR 
HRRO( IMA IN ) = SNGL( 3 Sll“ R ) 

CALCULATE IMAGINARY ®V»T OF HR ( W , T WflP I F > 

I FUN= 2 

FORM THE SEQUENCE OF PARTIAL INVERSE HANKEL TRANSFORMS BY INTEGRATING 
BETWEEN LAMBDA SUB(I) AND LA M3 DA SUft(I»l)f FOR r = I TO NARRAY (THESE 
LIMITS ARE IN THE ARRAY DHIfJTL). 

INITIALISE the VARIABLE USED TO ACCUMULATE the results of the 
NUMERICAL INTEGRATION OVER EACH L A M 3 0 A SUB-RANGE* 

dint r -n . no 
DU 3 1 = 1, NARRAY 

FORM TH F LOWER AND UPPER INTEGRATION LIMITS 
DA=DHINTL( I ) 

DB=DHINTL( I ♦ l > 

INTEGRATE (THESE VALUES ARE PASSED TO THE FUNCTION DF C T THRU COMMON). 
RESULT OF THE INTEGRATION IS STORED IN DA3RAV( I ) FOR FUTURE USE IN 
The EULER CONVERGENCE ACCELERATION ALGORITHM DTEUL. 

CALL ■JGAUS5(DA,DB,r c UN,DRER,DAP.SER,DARRAV(I),DERR t IE») 

STORE ERROR Flag 

I ERFLG( 2, I MAIN, I ) = r ER 

NOW SUM The INTEGRALS FOR COMPARISION WITH THE ACCELERATED SUM TO 
BE CALCULATED IN SUBROUTINE OTfUL. 

) DINTI=DINTI*DARRAY(I) 

NOW USE The EULER CONVERGENCE ACCELERATION ROUTINE TO BEST ESTIMATE 
THE REAL PART OF THE TOTAL INVERSE HANKEL TRANSFORM AT THE 
FREQUENCY TUOPIF. 

CALL 0TEUL(QFCT,CSUMI,NAR3AY,1.E-I4»IERI»NI) 

Oh I R 0 ( IMAIN) = 0SUMI 
HIR0( IMAIM)=SNGL( DS'JMt) 

F = 0 M E G A ( I MA in ) / S. ,?B ) 1 33 

BOO WR I T E(o, NO 10 ) F , DHR R 0 ( tMAIN),IERB,NP ? OHIPil(lMAlN),IERI»NI,OINTA, 

~ D I N T r 

90 10 F0.RMar(Y<4,F7.O,TIw,2(D23.1b,2X,Il,2X,I3,l2X) f /,TIO,*SUM = *,023.lG t 
STS3, • >UM= » , 02 3 . I 6 , / ) 

FN3 OF DEVELOPMENT STAGE CODE 

BLOCK 3 - FORM HR(R,N— DELTAU) NEFOED FOR IDF T. USE SPLINE FUNCTION 

INTERPOLATION ON THE RESULTS OF BLOCK 2. 

NSPLr 77 

C FORM SPLINE COEFFICIENTS FOR THE REAL PART OF HR(A,0) 

CALL SPL I N r ( NS PL , DM EGA, HARO, HR SB ,HR SC * HR SO ) 
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ORIGINAL PASE'IS' 
OF POOR QUALfTY 


C rnin t,PL!N‘ f f f IC t »■ N TS FO^ THf IMAGINARY PART OF HP(P,0) 
c all -."L INS ( NSPL ,UMFGA «1,H I SO ,H ISC »HI SO) 
c r IJP >i*MC YEPS E(|p use IN I D f r calculations 
r s A M »• - " . F - s 
N 0 F T - n 3 4 
TVIUP I = 2.* 3 . 14 ISO 3 
U:j = TWUP T / < r SA**p$NnF T ) 

I STOP = N[)F T/2, 1 

C T M ‘ SPrCHUH OF THF ([JST (N/2*l) FREQUENCIES 

DO 4" ! = i. I STOP 

3aBF«') = WOFFLOAT( 1-1 ) 

SO C F«a ( I ) = C MPLX ( SE V ALP ( N'.PL , * AfiFPO ,OMEr,A .HP AO ,np S3 ,MP SC ,HP SO) , $E V AL I 
*< NSPL. A A [)F-i. OMF 0A,HI :*n,HtS‘l,K ISC t MI SO) ) 

C NOW FOAM THF CONJUGATE PAST OF THF SPFCTPUM 
rSTA»T=ISTJP*l 
DO SO I=ISTAPT ,NDFT 
SO CFHR ( I ) = CONJG( CFMP(N0FT *2-1) > 

c 

C (LOCK 4 - CALCULATE Th? INVERSE FOURIER T A AN S F 0 < M OF THE MAGNETIC 

C INTENSITY. USES AN I OF T TO SIMPLIFY THE NUMERICAL 

C PA OC E DUP E • 

c 

n£LTA=l./( TSAMPOFLOA Y(N0F T ) ) 

KO DF = - 1 

call FFT (KODE ,NDF T,DFLTA,CFHP , ASO) 

C “PINT P r: SUL T S - TIME DOMAIN MAGNETIC INTENSITY 
WRITE (4,2000) 

3000 F,1PMAT(//* TINEIhILLISEC ) INDUCTION (TESLA)*) 

DU TO 1 = 1, 2CI 
TIMCzr.ooSSFLOATC t-1 ) 
cfh-u I ) = c FHP( i )& c*un 
to UR I TE(S, 2010) TIME,CFH>.( I ) 

2010 F0PMAT(4X,Fa.I,r22,F8.S»T<tl,Ell.4) 

DO SO I = 2 2*, 100 1,5 
TIME = .COS«=FLOAT( 1-1 ) 

C F HP ( I ) =C FHP ( I )$CMUO 
SO WAIT F(s, 2010) TIME,CFHP(I) 
c ‘’PINT CPPOP FLAGS 
W.PI T E ( s , 4 SO 0 ) 

4500 FCJPMAT( * 1EPR0P FLAGS FPOM PEAL INTEGRATION:*/) 

no 2 is i=i,napray 

2 14 WPI re (.,,4 00 1) (IE pFlG(1*II, I), 11=1,77) 

4 501 FOP MAT ( 2X , 77 ( II) ) 

WP I T>- (4,4502) 

4502 FOPMATC 1EPR0P FLAGS FROM IMAG INTEGRATION:*/) 

DO 217 I = l , NAP RAY 

217 UP I T F ( S , <» 50 l ) ( IEPFLG(2,II,I),II=l,77) 

C 40 CONTINUE 
si STOP 
END 

REAL FUNCTION S6VALP(N,U,X ,y,b,c ,D) 

C COPIED SUHPPOGPAM SEVAL FOLLOWS 

PEAL U,X(N),Y(N),3(N),C(N),D(N> 

DATA 1/1/ 

IF(I.GE.N) 1=1 
IF ( U.LT. X ( I ) ) GOTO 10 
IF(U.LE.X( 1*1) ) GOTO 30 
10 1=1 
J = N* l 

20 Ks(I»J)/2 

IF(U.LT.X(K) ) J = tC 
IF (U.UF . X ( X ) ) I = K 
IF(J.GT. I ♦ 1 > GOTO 20 
10 DX=U-X(I) 

SEVAL R = Y(I)*DX$(a(I)*DX*(C(I)*DX«=D(I>>) 

PE TURN 
END 

PEAL FUNCTION SE VALl (N,U,X ,Y, 0,C ,H) 

C COPIED SUHPPOGPAM SEVAL FOLLOWS 

PEAL U,X(N),Y(N),H(N),C(N),D(N> 

DATA 1/1/ 

IF(I.GE.N) 1=1 
1F(U.LT.X(D) GOTO 10 
I F ( U. LE • X ( 1*1) ) GO r 0 30 

10 1 = 1 
J = N* l 

20 <=( I*J)/2 

IF ( U.LT. X (K ) > J = K 
IF(U.GE.X(<) ) I=K 
I F ( J.GT. 1*1 ) GOTO 20 
SO OX=U-X(t) 

SEVALI=Y( I )*DX~(B( I ) * DX-(C ( l ) *DXOD( I ) ) ) 

RETURN 

END 
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c 


c 

c 

c 

c 

c 


subrouti SPL tNE(N,x,y,fl,c»n> 

CU^'lEO SUBROUTINE SPLINE F OL LOW 3 

t -M T r ’ -> N 

4 E 1 L X ( N > , Y ( N > , H ( N ) , C ( N ) , 0 ( N ) 

[ N T f f P N M l • IB, I 

»C4L r 

NMl=N-l 

IF ( N.LT . 2 ) 3C TURN 
IF(n.lT. 3 » ' i 0 TO 50 
D(i)=x(2)-x(i) 

C(?)=(V(?)-V(1>)/0(1) 

no ioi = a , nhi 

um = x< i «- 1 > -x ( i ) 
i) < i ) = j.«(o(i-i)*o(in 
c< > = (y{ i*i)-y{ t > > ✓ ot r > 
c < n = c< r ♦ i > — c < t ) 

10 CONTINUE 
b< i >=-n< i) 
h( N ) r-n( N-l ) 

c(i)=o. 

C(N)=0. 

IF{N.fQ.J) 50 TO r 

C< 1 )=C( })/( X(<*)-X(2) ) - C( 2>/(X( 3 » - X ( l) ) 

C(N)=C(N-1 )/<X(N)-X(N-2)) - C(N-2»/( X ( N- 1 ) -X ( N- 3 ) ) 
C( I )sC( n*tl( 1 )«*2/(X(«,)-X( 1 ) > 

C (N > =-C ( N ) *D(N-1 )»i=2/( X(N)-X (N-l) ) 

1 5 00 20 I = 2, N 

Tr 7{r-l>/BCI— 1) 

RU) = MI>-T«D(t-ll 
C ( I ) -C < l)-T*C( 1-1 ) 
ao c o n t four 

C(N)-C{N)/R(N) 

00 30 I'l= l, NHI 

r-N-rri 

c ( n = <c 1 1 )-o< i >*c ( i *i ) )/»( i ) 

30 CONTINUE 


»<N) I (Y(NI-Y(MNim/n(N*U) 

00 4 0 I = 1, NHI 

b c i ) = (y(ih)-v(i n/o( : ) 
0(1) = (C(I*1)-C( t n/n(t) 
C( I ) = 3 • « C ( i ) 

CONTI MU 1 -: 

C ( N ) = 3 .*C(M) 

0(N)=D(N-1> 

RETURN 

n(I) = (Y(2)-Y(l))/(X(2)-X(U) 

C( 1)=0. 

0( 1 )=0. 
f*(2)-0(l) 
c ( a ) = o . 
o(a)=o. 

RE TURN 


n(NHl )*(C (NHI ) » 2.«C(NI) 

cm*(c(i*i)* 2 .*:(i)) 


END 

SUBROUTINE FFT(K00F»N»DELT4»X*R) 

file: fft fortran 

PROGRAM Tn IMPLEMENT THF FAST FOURIER TRANSFORM UH F N THE NUHBER OF 
DATA POINTS IS AN IMT'C.cr p OWFR OF TWO. TAKEN FROM PAGES 2S4-2S3 
OF "METHODS OF DISCRFTE SIGNAL AND SYSTEM ANALYSIS”, BY JONG. 


COMPLEX x<n>,u,xi,cnplx 
ia =c 

M 1 =N 

s Na=Ni/a 

IF(Na : .' 2 .NE.Nl> GOTO 100 

I a = I i » 1 

N l = N a 

IF(Nl.r,T. DGOTO 5 
”N = S. .’ 3 .U 3 S/N 

l = n / a 
t»i=ia-i 
K I - 0 

00 30 I S = l » I R 

is oo an i - 1 ,l 

K»K I * l 
KPL=K*L 

AMzKD l TR ( Kl/a«®IRl, I P ) 
IF(AH.NE.O.) GOTO 13 
X 1 = X (KPL > 

GOTO LO 
13 ARG = AHC=PN 
C = C 0 S ( A? G ) 
S=-KOOFOSIN(ARG) 

W-Chplx ( C , S ) 

X l - W * X ( X a L ) 
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I I X(KRL> = «<0-X1 
X(*);X(<)»X1 

,’0 V l - < t * t 

X- i :K 1 »l. 

IF(ki.Lt.N) r.oT'i 15 

v. i * o 

r » i = r •' l - i 

!'■> l-l/q 

00 K:1,N 

K 1 CK ">1 T UX-1 , l R ) »l 

iMKi.ir.x) i,n r n «.o 
X 1 = X ( K ) 

X(K) r*(H 1 ) 

X(*l )=X1 
LO CONTINUE 

t p ( DELTA . E Q. 1. ) RETURN 
no 50 K = 1*N 
>o x(k ) = 0 ?t.rasx(K) 

’ E TURN 

100 UP l TE(*>, 101 ) N 

101 F(JRMAT(//,* S5S M= •»!<»»• 15 MtlT a POUEP OF 2, FFT RUN ABORTED 

I •- 1 

•0 E TURN l 
END 

FUNCTION K3tT«(K,IP> 

KBITR-0 
X 1 -K 

0(1 1 [ - 1 1 1 P 
K2=Kl/2 

KBITRs2®K5ITR*Kl-2®K2 
1 *1=K2 

RETURN 
END 

DOUBLE PRECISION FUNCTION DFUN(DX) 

C SUBROUTINE TO FORM THE INTEGRAND IN GAUSSQ. SINCE GAUSSQ IS USED FOR 
C NUMERICAL QUADRATURE WITH TUO DIFFERENT INTEGRANDS, WHICH INTEGRAND 
C IS EVALUATED BY FUN IS DETERMINED BY THE FLAG IFUN, PASSED TO FUN 

C A V THE MAIN PROGRAM THRU COMMON STORE. 

IMPLICIT COMPLEXES (C), REAL®B (D) 

C0“plEXS 1E, CO/<.aiO-3,O.QO)/',CQNE/(l.DO,O.DO>/,CTUO/(2.DO,O.DO>/, 
■50CMPLX 

COMPL-X curnt 

reals 1 Dmijsid/ 21.9G5‘»B00 / , OH D3/1.5RD — 3/ »D5/2.TaD— 3 / ,DCOMP( 2), 

S n A R R A Y I 1 0 0 ) 

‘-QUIVALFNCE ( CHR , DC0MO( l ) ) 

COMMON UAPRAY, CURNT, IFUN, R,TUOPIF 

C Ff)RM COM*»L r X® l £> VERSION OF THE COMPLEX®a CURRENT PASSED TO IFUN THRU 
C THE COMMON VARIABLE CURNT 

CURNT!)=UCMPLX(D3LE(REAL(CURNT>),D3LE(AIMAG(CURNT))) 

I F ( D X • N 6 • 0. DO ) GOTO 10 
C LAMftDArO; I N T F GR AN 0 = 0 
DEUNsO.DO 
RETURN 

C EVALUATE THE NUMERATOR (WITHOUT THE KERNEL) 

10 DTi = ;ifXP(-DX®0HD3) 

DT 2 =((nri*l.DO)®Ori*l.DO)®OEXP{-DX®DG>®DKPRLK(DX) 

C IF FR £ QUE NC V =0 , ALMOST DONE 

I F ( T WOP I F , NE • 0.) GOTO 20 
C FPEQUENCY is zero 

CHR=nCMPLX(0r2/2.DO,0.00l®CURNTD 
GOTO 10 

C neither lambda nor FREQUENCY is jfro 

RO CSQ=CDS3RT(DCMPLX(DX®nX,DBLE<TW0RIF )®D*USIG) ) 

CT=COEXP(CD®CSQ) 

C T I = C ONE /C T 

CCOSM=(CT*CTI)/CTWO 

CSINH=CCOSH-CTI 

CR AT IO=OCMPLX( DX ,0.00) /CSQ 

C D E NO M=C ONE ♦ ( CCQSH.CR AT 10®CSINH)/(CSINH/CRATI0*CC0SH) 

CHR=DCMPLX( DT2,0.30)®CURNTO/CDFNOM 

C NOW CALCULATE integrand 
30 DRTL=nX*DBLE<P ) 

D Jl = D J IL K ( Dp TL ) 

1 F ( IFUN . F Q . I) DFUN=DCn“*’( 1)®DJ1®DX 

I F ( IFUN . E Q • 2) 0FUN=nCClMP(2)SDJlSDX 

RETURN 

FNO 
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douri.: ^FCiMnN function dkprlk ( dlmbda > 

C SUBPROGRAM T'l EVALUATE Tuf FUNCTION < PH I M fc ( L a M R 0 a ) . THC ARGUMENT 

c pl mji r> a must Mt- pfflL- 0 .. u s 6 ■> the result that the ant r ofr i va r ive 

c n= j i ( x > * x is e jijal rn 

C P t * X * { J1 ( X ) * MO( X ) - Hl( X ) * Jimi / 2 

c Ci H A?..,u“t nts dLMitnaopx .le. 0 . the difference is evaluated 
c directly. for clm'idasrx .or. a, some s impl tf rc at r om rs 

C UfILI’ l 2 (SEE the PAPEJ DATED JULY 23). 

C » 5 ;f: "HIGHER T R AN SC -E N 0 E M T AL F UNC T J ON S” , VOL. 2, P. 39, 

C EQUATION (5) WITH NU= 1 

C PEgUTPES E U NC 1 1 fl N SUE.PR OC.tr AMS flHGL K ( D X ) » OH 1L K ( 0 X ) , 0 AO SM ( DK ) , 

C DAXSM(llX), AND THE REAL-0 IMSL SURxUUTINE “MftSJN. 

C M, 10 1 F I f. 3 SEPT 19 AS AT USU TO REPLACE MMBSJN WITH MY E.P. FUNCTIONS 
C OJOLK AND D J ILK . 

C 

IMPLICIT P E AL* A (0) 

PEAL- 3 IH* 1/ 3. 14 l SO ESS 35.3979/ , DR 1 / . 0 0 3 1 7 5 DO/ ,DR2/.0254DD/,DJU(2), 
*!IJL ( 2 ) 

DATA IFLAr./-l/ 

C C H c c K TO S' F IE ROUTINE HAS SEEN CALLED BEFORE 
IE (I FLAG . NF . -1) StlTO 10 

C NO PREVIOUS CALL - PE9EU»« INITIALISATION 
DC ON ST rS. PO* OP I/( DPS- DR 1 ) 

C RESET FLAG if) INITIALIZATION WON*T HE PERFORMED AGAIN 
IFL«:.--1 
C 

C CALCULATE V PM IMF ( LAMBDA) 

c 

10 DK PRLX- ■' . n ) 

C C H f C < Til St ‘ IF ARGUMENT = 0 

I F ( [I L M R A . ?Q. 0 .DO) RETURN 

C RESIN RV E Or M I N(j IJPPEP AND LOWER INTEGRATION LIMITS 
DU = DRS*[)Lm:iDA 
DL = Di> 1*DLM IDA 

C FORM HFSSEL FUNCTIONS OF THE FIRST KIND AT THE UPPER AND LOWER 
C INTEGRATION LIMITS 
l)JU( I ) = DJOLK(DU) 

/)JU(2) = PJILK(DU) 

!) JL ( l ) = P JOLK ( DL ) 

0JL(2) = ::jlLK(0L> 

C FORM DIFFERENCE AT THE UPPER INTEGRATION LIMIT 

C DECIDE UHTCH TECHNIQUE TO USE Tf’ EVALUATE J l ( DU )*H) ( nu > -H l ( OU )* J 0 ( DU ) 
IE (DU .GT. 0.00) SOT.) 20 
;)IFFU = DJU(2 >*DHOLK< DU ) — DH ILK ( DU ) * D JU ( 1 ) 

GOTO 30 

.?•.) DIEFU = DJU( 2 )* DAOSM( DU)-nAiSM( DU)*DJU( 1 ) ♦ 2 • DO / ( DPI* DU ) 

c eqjm DIFFERENCE AT THE: LOWER INTEGRATION LIMIT 

C n F c IDE WHICH TECHNIQUE TO USE to EVALUATE Jl( DL )*H7( DL)-Hl( OL)*JO( DL ) 
10 IF ( DL -GT. 3. DO) GOTO t»0 

DIEFL = DJL(S )* PhOLK ( DL)-0MILK( UL >*.DJL( 1 ) 

GOTO SO 

40 niFFL = 0JL(2)« : nA0SM(nL)-DAlSM( DL)-DJL( 1) ♦ 2 . DO / ( 0 P I * DL ) 

C EVALUATE KPRIME 

SO DKPRLK* DC ON ST* ( 0R2-DI FEU-DU* DIFFL )/DLM3DA 

return 

ENO 

FUNCTION DHOLK(DX) 

C n 0J RL E PRECISION SURPRQGRAM TO EVALUATE THE STRUVE FUNCTION OF 
C ORDER ZERO FOR DOUBLE PRECISION ARGUMENTS DX .LE. 0 USING A 
C SUMMATION OF weighted chehysmfv polynomials. 

C REQUIRES The SUBROUTINE DCNPS (NASA LEWIS COMPUTER LIBRARY) TO 
C FORM AND SUM THE polynomials. 

c rfe: ”The s oc cial functions and their approx imations", 

C VOL. 2, P. 370, BY YUDFLL LUKE. 

C 

IMPLICIT realms (0) 

REAL-0 B(l>)/1. 00215043009012, -1.03969292031309,1. 50 2 369396 10293, 

*-. 72405115 302121 0,. 109 55 327371O931»-.O3O07O52O22 900, 

* 3 37. ->*,14 4 7 3 7SI 94 0 -5,-26. 945014 3 1260 2D- 5*1. 637 46169261230-5, 
*-7024. 4400 S 002540— 10* 30 2. I 5 9 3 1 0 3 l 5 3 0 - 1 0 » - 9 • 6 3 2 6 6 4 4 9 5 D— 1 0 , 

*. 257 » 3370 3 9D- 10, -5 3 3. 5374 'll) - 15*11. 50332 D- 15/ 

C CHECK FlJR ARGUMENT OUT OF RANGE 
IE{0AR$(0X) . G 7 • 0.00) RETURN 

C FORM The ARGUMENT FOR THE CHEBYSHEV POLYNOMIALS 
OXO9=DX/0.DO 

DARG=2. DO *0X00*0X03 -1.00 
C SUM THE POLYNOMIALS 

CALL flCNPS( UHOLK , DARG ,3 f IS ,09) 

C COMPLETE The evaluation AND RETURN 
DHOLK =l)HOLK* DX 00 
»9 RETURN 
END 
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sui'-iijur i fit-; nCNPS(Y,x,c,N,*) 

c 

;; I » r: | C(N) 

double ” fi c I s i on c.X|X,hO(M1,h.’, a*g 
c 

C TEST l)F DIMENSION 

if(oahs(x> .r,T. i .o?)re turni 

f F ( N >1 , t , 2 

1 RETURN 
C 

2 IF(N-2> 1,4,4 
1 Y = C ( 1 ) 

R E T U ' N 

c 

c initialization 

<* ARG = K*X 

hi = 0. an 
no so. oo 

c 

do ■> r = i,M 
k = n-i 

H2 =H1 
H1=H0 

- H0=flRG*Hl-H2»C (K+l > 

Y=0.500*(C< 1)-H2*H0) 

RETURN 

END , 

FUNCTION UHIL'K(DX) I 

C DOUBLE PRECISION SUBPROGRAM rj EVALUATE THE STRUVE FUNCTION OF 
C OPOfcP ONE FOP DOUBLE PRECISION ARGUMENTS DX .LE. ft USING A 

C SUMMATION WEIGHTED CHE BY SHE V POLYNOMIALS. 

C REQUIRES Th= SUBROUTINE DCNPS (NASA LEWIS COMPUTER LIBRARY) TO 
C FORM AND SUM THE POLYNOMIALS. 

C p F : "THE SPECIAL FUNCTIONS AND THEIR A P PR OX I M A T I ON S* , 

C VOL. 2, P- 170, BY YUDELL LUXE . 

C 

IMPLICIT REALSft (D) 

RE AL* ft C(19)/.557Aft91 4444814,-. 11 1 ft A 32 5 7 24549 ft , - . 1 6 1 37 9Sft 1 2 52 009 , 
*.322>491207240S9,-.l45A1432lb72442,3292.4TT 3993740-5, 

*-4oO. 17 21 4 209 3573 0-5,44. 3 4 7 0 4 1 ft 3 J 1 3940-5 , - 3 . I 4 20 995293411 7 D-5, 

*. 1712371993 30 035 0-5, -741. 4 9ft 7 005204 D- 10, 24. 13374707050- 10, 
*-.7435519 3 *50-10, 1904.70414D-15»-40.522910-15/ 

C CHECK FUR ARGUMENT OUT OF RANGE 
I F ( 0 A 55 ( D X ) -GT. ft. 00) RETURN 
C FORM THE ARGUMENT FOR THE CHE 5Y SH E V POLYNOMIALS 
DXOft-OX/ft.OO 

0 ARG= 2. 00*0X05*0X03-1 . DO 
C SUM THE PULYNOMIALS 

CALL OCNP S( 0H1LK , OARG.C , 15,99) 

T9 RETURN 
END 

FUNCTION OAOSM(OX) 

C DOUBLE PRECISION SUBPROGRAM TO EVALUATE THE SUMMATION ASSOCIATED 
C WITH THE STRUVE FUNCTION OF ORDER ZERO FOR DOUBLE PRECISION ARGU- 

C MFNTS ax • G F • ft USING A SUMMATION OF WEIGHTED CHEBYSHEV POLYNOMIALS. 

C REQUIRES The SUBROUTINE OCNPS (NASA LEWIS COMPUTER LIBRARY) TO 
C FORM AND SUM THE POLYNOMIALS. 

C P‘-F; '«THf 5PECIAL FUNCTIONS AND THEIR APPROXIMATIONS", 

C V'JL. 2, P. 171, BY YUOELLLUKE. 

C 

IMPLICIT REAL* ft (0) 

REAL* 5 0(20)/. 992 517275744239,-4 94. 8912ft l 13B4250-5, 

* 13. 20510 37 ft TO 3 71 0-5, -1.04 1258 25 2544 140-5, 93 19. 32942 845250-10, 

*-122 1.0445444977 0-10, 1 39. 40ft 111 l ft D- 1 0 , - 34 . 4 3 5 3 22 5404 D- 1 0 , 

*7. 1 l 19 10 171 ID- 10, -l.b2ft ft 7441 37 0-10, .40454807230-10, 

*-. 13 r ' I SO 47940-10, 5120. 0524 ID-15,- 94 2. 0207 0-15, 298.47 947 0-15, 

*-93.7 2 414 0-15, 3 I. 9 371 2 0-15, -12. 07930-15, 4. 439210-15,-1.47 85 9D-15/ 

R E AL* ft DP I 02/1 .5 70794 3247949/ 

C CH'CX FOP ARGUMENT OUT OF RANGE 
IF(OAftS(OX) .LT. ft. DC) RETURN 
c cfjJM THE ARGUMENT FOR ThE CHKiYSMPv POLYNOMIALS 
08UX-ft. OO/DX 

DAR G = 2. 00 *030X* D8DX-l.no 
C SUM the POLYNOMIALS 

call dcnps( daosm, ua->g, 0,20 ,99) 

C COMPLETE the evaluation and return 

DAOSMr DAOS m/( op I 02* ox ) 

49 RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 

c 

c 


tunc r I rjN on ism( ox ) 

double precision subprogram tb evaluate the summation associated 

with the STRUVE FUNCTION of ORDER ONE FOR DOUBLE PRECISION APGU- 
mfnts*'cx .ge. .a using a summation of weighted chebyshev polynomials. 

REQUIRES THC SUBROUTINE DCNPS (NASA LEWIS COMPUTER LI3RARY) TO 
FORI AND SUM THE POLYNOMIALS. 

R<:f: "THE SPECIAL FUNCTIONS and ThEIK APPROXIMATIONS"* 

VOL. 2, P. 3 7 1, BY YUDELL LUXE. 

IMPLICIT PEAL- A Cni 

REAL-* E ( 17 >/l .00 TG7647 29 3*66 ,7*0. 31605124*2570-5, 

S- 7. 9 4 393 32645 l 905 D-5. . 26629 5 3 9 3 3 8 2266 D-5 , - l ft *4 . I 1 3 7 7S 3405 D- 10 , 

2 194.901495*3940-10, “26. l26l9ft990SD-10,4.2362690l04D— 10, 

TO OS 1555 Un-10,. 167997 30 OSD- 10,- 3007 . 193210-15,985.43090-15, 
R-266. 357940-15,76. 4 50 35 n- 13, -23. 129610-15, 7. 3 3212 D- 15, 

S-2. 423340-15/ 

REAL-* 0PIO2/1. 5707063267949/ 

CHECK FOP ARGUMENT OUT OF RANGE 
I F ( 0 ADS ( D X > .LT. ft. DO) RETURN 
FORM THE ARGUMENT FOR THE CHEBYSHEV polynomials 

0 3 D X = 3 » 0 0 / 0 X 

0APG=2. 00-03DXO030X-1 .DO 
SUM THE POLYNOMIALS 

CALL f)CNPS( 0A1SM, 0ARG,E,17,90) 

COMPLETE Tut EVALUATION AND RETURN 
DA 1 M = D A 1 SM/OP 102 
TO RETURN 
END 

DOUBLE PRECISION FUNCTION DJILK(DX) 

DOUBLE PRECISION SUBPROGRAM TO EVALUATE THE BESSEL FUNCTION OF THE 
FIRST KIND OF O-tOFR 1 AT THE DOUBLE PRECISION ARGUMENT DX. FOR 
ARGUMENTS DABS(DX) .LE. 8 A TRUNCATED JACOBI SERIES OF CHEBYSHEV 
POLYNOMIALS IS USED (SEE THE FIRST REFERENCE BELOW). for ARGUMENTS 
DABS(DX) .GT. * A TRUNCATED ASYMPTOTIC SERIES IS USED (SEE SECOND 
REFERENCE below). 

RFyuiRfS THE SUBROUTINE DCNPS (NASA LEUTS COMPUTER LIBRARY) TO FORM 
AND SU“ THE CHEBYSHEV POLYNOMIALS. 

REE: "THE SPECIAL FUNCTIONS AND ThEIR APPROXIMATIONS", 

VOL- 2, P. 332, BY YUDELL LUKE. 

>ff: "CHEBYSHEV SEPIES FOR MATHEMATICAL FUNCTIONS", NAT. PHYS. 

LAB. MATH. TABLES, VOL V, P. 33, BY C.W. CLENSHAW. 

IMPLICIT REAL-3 (D) 

REAL' 1 * D B (15)/. 64335 *77060 5264 00,-1. 1 9 1 30 1 1 b 0 54 1 2 2 DO , 

* l .2*7 994 0 9**5 76* DO, -.66 144 X4 1 3 454 i'OO, . 1 7 7 7091 i 72 3972300 , 

P -29 17. 5^24*061542 0- 3, 3 2 4. 0?7D1 ft 2 633*60-5, -26. 04443 *9 34 35 310-5, 

* 1 • 3**70 1 92 39932 ID- • ,-76 17. 537* '5400 3D- 1 C, 294. 9 70 7 007 27 3D- 10, 

! S -9. 4, 42 129316U-1 0,2 5 231.23^640-1 5, -577. 74 04 20-15, 

911. 3*5720-15/ 

REAL** DP 1 ( 15 >/2.00l *260*1 7200 JDO » ft 9. ft 9ft 9 * J 3035)41 D-5, 

> - 39* T2. 34 3004389 i;j- 10, 5 I 7. 76 3 396 06440-10, - l ft. 71 *9074 911 0- 10, 

*331 b*.9ft h 6U0-l 5, -5704. *o 364D- 13, 469. 9 19550- 15 ,-46. 342240-1 5, 

- 5 . 4 5 267D - l 5 , - 7 2 2 1 2 . 0 — 20 , 1066* . D— 20 , — 1 7 3 l . 0 - 20 , 30 5 . D— 20 , 

3-53.0-20/ 

J E«LR3 D01( 15)/ 1 355.55741 3907070-5,-9.6277235491570*0-5, , 

-91 33. 61^26795350-10, - 209. 5R731 3 *40 *0-10, ft. 229193 3 277 D- 10, 

--46*63. S36**D- 16, 3615. 2187 9 0-13, -526. 431570-15,35. 96 7770-15, 

4=- 4. 56 12 511-15,6 50*j.D-20,-l 926 9.0-20,17 68.0-20,-3 23.0-20, 

*65.0-20/ 

DATA DP 1/3. 1415926535 ft 979/, IFL AG /-!/ 

CHECK TO SEE IF ROUTINE has been CALLED before 
r F ( I FL AG .NE. - 1 ) GOTO 5 
PERFORM INITIALIZATION AND RESET FLAG 

ocoNsr=nsqRr(2.Qo/oPi) 

D JPI 04= .3. DO- OP 1/4. DO 
DPI ( 1 ) = DPI( l )/2. 1)0 
031(l)=DOl(t)/2.00 
iflag=i 

3 OJILK=O.DO 

CHECK FOR ARGUMENT cUUAL TO ZERO 

1 F ( 0 X .r ; q. 0.00) RFTURN 

DECIDE IF JACOBI OR' ASYMPTOTIC SERIES WILL RE USED 
I F ( riABSI DX ) .GT. ft. DO) GOTO 19 
USE JACOBI SE 3 I F S 

FQRM THE ARGUMENT Fqr THE CHFBYSHEV POLYNOMIALS 
nXC8=DX/3.D0 

0 A R G = 2 • 00*0X08*0X0 ft- 1 . 00 
SUM THE POLYNOMIALS 

CALL DCNP S( 0 J1 LK , OAP.G, OB, 1 3 ) 

FINISH THE CALCULATION AND RETURN 
OJ I LK 6 Ox 0 ft— 0 J 1 L K 
99 RETURN 

USE “QDIFIEn ASYMPTOTIC SERIES 
10 D30X=8.DQ/DX 

OAR G p 2. DO * D ft OX SO ft OX - 1 . DO 
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c 

c 

c 


c 

c 

c 

c 


GU" T Mf pULYNUmj OLS FOR Pl(X) 

CALL !»C'|P S( DP , dari,, UP 1 1 IS) 

5I)m r«h P(lLYNilM[4LS FOR Ul(() 

CALL I'C NPS( 93, BA5 U , 1)01 ' 1 5) 

o o - o j d n x 

Ffl'M TH< ASYMPTOTIC EXPRESSION 

DJ ILK-( DCI)NST9< OPillCOSt DX-D 3 P I 04 > - DQ9 OS I N( DX-0 JP 104) ) )/DSQR T( OX > 

RE TURN 

FMD 

DOUBLE PRECISION FUNCTION ojolk(dx) J 

DOUHLF PRECISION SUBPROGRAM Til £ V ALUA T E THE BESSEL FUNCTION OF THE 
FIRST KIND 0 F fl A 0 E A 0 AT THE DOUBLE PRECISION ARGUMENT DX. F OR 
ARGUMENTS DAHS(DX) .LT. A A TRUNCATED JACOBI SERIES OF CHEBYSHEV 
POLYNOMIALS Is U S f D (SEE the FIRST reference BELOW). FOR ARGUMENTS 


c 

C 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


0 A a s ( DX ) -r, E. a A TRUNCATED ASYMPTOTIC SERIES IS USED (SEE SECOND 
AND Third RFEfRfNCE BELOW) . 

REQUIRES THE SUBROUTINE OCNPS (NASA LEWIS COMPUTER LIBRARY) TO FORM 
AND SUM TH£ CHF !1Y SHF V POLYNOMIALS. 

RIF: "T'lF SPECIAL FUNCTIONS AND THtIA APPROXIMATIONS", 

VOL. 2, P. 151, BY YUDELL LUXE. 

pee: m CMF:!YSMEV SERIFS FOR MATHEMATICAL FUNCTIONS", NAT. PHYS. 

L Ail. MATH. TABLES, VOL V, P. 12, BY C.W. CLENSHAU. 


IMPLICIT REAL98 (D) 

REAL- a 0M1S)/. I 5 7727 4 71474 890 DO, -872. 3 44 23528522? D- 5, 

9. 26517 8 613 20 J334D0. -.370094*9 31 7 26 4900, .158 04 7102332097 DO, 

9- 3484. 37 694 H4084D-5 ,481. U 8006944 7600-5 ,-46.062 6 166 20 62 7 50-5 , 
9 3. 24 40 3288. MOO SID- 5, -17619. 4k 9,0 776215 0-10,740. 81635924 19-0 -10, 
9-2 6.7125 35 3056 0-10, 7 1 416. 9-, 3 140-15,-194 3. d 3469 D- 15, 

941.24 3210— 15, — 75185.0 — 20,12 22. D-2 0 , - 1 7 . D-20/ , 

9 DPO( 14) / 1. 9 9 8920491 69 50 4 00 , -5 3. 6522 04611 3212 D- 5, 

9 30 751 .84 7175 1947!)- 10, -5 17. 0 5945 J 7 6060- 10, 16. 30 64646 35 2 D- 1 0 , 
9-786 40. 91 377 D- 15, 5141. 262 3 9 0—1 5, -430. 45789D-15, 4 3. 26 5960- 15, 

9- 5.0690 50- l 5, 6 74 3 1.0-20,-100 l 2. 0-20 ,16 31. D-20, -288. 0-20/, 

9DQ0< 14)/-. Ull. 17092 1047400-5,6.8 38519942611650-5, 
9-74l4.4'»84llOb060-10, 1 7 9. 7 245 724 797 D-l 0,-7. 27 19159 3690- 10 , 

9 42 20 1- 2 1 9040-1 5, -320b. 7 47 42 D- 15, 300. 61 4 51 0-1 5, -3 3. 363280-15, 
94. 254230-15,-60999.0-20 , 96 b2 . 0-2 0 , - 1 66 9 . D-2 0 , 311.0-20/ 

DATA DPI/ 3.14159245 3 4 8979/ , I FLAG/— 1/ 

CHECK Til Sf: f IF ROUTINE HAS BEEN CALLED BEFORE 
I F ( I F L A i j -NC. -1) GOTO 5 
PERFORM INITIALIZATION ANO RESET flag 
DC 0NST = I)SORT(2.OO/0PI> 
nPID4=OPI/4.00 
0PC( 1 >=0P0( l)/2.00 

noo( 1 )=DOO( D/ 2.00 
IFLAli-1 
5 I) J OL x = 1 . OO 

check for argument equal to z c Rn 

I F ( 0 X . r'J . 0.00) RETURN 

DECIDE IF JACOBI OR ASYMPTOTIC SERIES WILL BE USED 
IF(0A.<1S(0X) .gt. S.DO) GOTO 10 
USE JACOBI SERIES 

FORM THL ARGUMENT FOR THE CHE <3Y SHE V POLYNOMIALS 
DXD8 = DX/a.t)0 

OA?G = 2 .00 9fJX0 8* 0X08-1. 00 
SUM THE POLYN'tMIALS 

CALL OCNPSI OJOLK , DARC.,08,1 8) 

4 9 .7 £ T U R N 

US C MODIFIED asymptotic serifs 

10 08CX=3.0G/PX 

OARG=2.009080X9080X-1.00 
SUM THF polynomials EOR PO(X) 

CALL DCNPS( OP, DAPG»DPO» 14) 

SUM THE POLYNOMIALS FOR QO(X) 

COLL OCNPS( 9Q»QA a G»nOO, 14) 

00 = 0 ,’*D? DX 


NOW FORM ASYMPTOTIC A p “ R OX I M A T I ON TO JO(X) 

njCLK = (DCflNST9(0P9 0COS(nX-DPID4) — OQ90SIN(OX — DPI04)))/0SQRT(DX) 


R ETU I N 
END 

DOUBLE PRFCISION FUNCTION DEC T( l ) 

FUNCTION TO EVALUATE TH*- INDIVIDUAL T^j^S IN THE SERIES OF MAGNETIC 
FIELD INTENSITIES ( E I r '< F H R ' 1 L i)F IMAGINARY) ARISING FROM THE 
INVERSE HANKEL INTEGRATION JETW'-f-N THE Ll“ITS OF LAMBDA SUB(I) 

AND LAMBDA SUM 1*1). SE^ THE P A P F R DATED JULY 25. 


IMPLICIT H E A L * 8 (D) 

R E a l 9 a darrayi nm 
common n a p p a y 

DFCT = :)A3R0r( I ) 

RETURN 

END 
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ORIGINAL RAGE F§ 
DC poor quality 


subr ji i r t n= rr r ul ( c c r • vii| m ax , eps , I r-R ,n ) 

C COPIED r‘UL ; A CONVERGENCE ACCELERATION SUBROUTINE 
D t “ c NS I ON r ( 1 >, ) 

DOUBLE PRECISION FC T , Su“ , Y , AM\ , AMP 
I F ( » A < U , 1 , £ 

1 IER=-l 

goto 12 

c initialize euler transformation 

2 icr = i 

i = i 

M 3 l 
N= l 

Y (I) =FC T (N) 

SU.M = y( 1 >*.500 
C \TAST EULER LOUP 
.1 J=0 

4 I = I ♦- 1 

IF < I-MAX ) 1 , 5, 1 2 
■5 N=I 

amn = fct( I ) 

DO S K = l,M 

AMP = ( AMN»Y(K ) )*.5D0 

Y ( K ) - AMN 

5 Amm = aMP 

C CHECK E OLE R TR ANSFORMAT ION 

IF(0ABS(AMN)-DABS(Y(M)))7,9,9 

3 M = M ♦ ! 

Y { M ) 3 AM*! 

AMN=.500*A“N 

C (J?U4fE SUM 

9 SUM=SUM.4“N 

t F ( A B 5 ( SNGL (AMN))-EPS-ABS( SNGL( SUM))>10,10,3 
C TEST END UF PROCEDURE 

10 J=J*I 

11 I E R = 0 

12 RETURN 

r N I) 

SUB-1 iUJT r NF (JGAlJSS(DLOWFR,r'UPPER»DFUN,DPELER»DADSER,DINT,DERR, 

* r f l a o ) 

c subroutine to estimate THE INTEGRAL of the FUNCTION DFUN between the 
c LIMITS OLOWC-t AND OWPPER. 

c PARAMETERS DWELT A AND D AD S£ 1 A* = THE USER REQUESTED RELATIVE AND 

C ABSOLUTE ERRORS IN the estimated VALUE OF THE INTEGRAL DINT. IF 

c OAbSERrf), (JNLY ThF RELATIVE ERROR CRITERION IS USED. COMPUTATION 
C STOPS WHEN EITHER THE ABSOLUTE ERROR CRITERION OR THE RELATIVE 
C ERROR CRITERION HAS BEEN ESTIMATED TO BE ACHIEVED, 
c PARAMETER uerr IS T h c RETURNED ESTIMATE of the absolute error in 
C DINT. PARAMETER I FLAG IS SET TO ZERO IF NO CONVERGENCE IS OBTAINED 

C OTHERWISE* I FLAG INDICATES the NUMBER of INTEGRAL ESTIMATES 

C CALCULATED MINUS I. 

C SUBPROGRAMS REQUIRED - D F UN , DGO U A D » C ONV T 

C REF; “GAUSSTAN QUADRATURE Fu»RULAS"i BY STROUD AND SECREST 
C (THE FUNCTION ARGUMENTS AND THEIR WEIGHTS WERE TAKEN FROM 

C THIS BOOK.) 

c 

implicit realms ( n ) 

R E AL*3 DGX 12(6 ) , 0GW12C6) , DGX l 6( 3 > , 0GW1 6( 3 > » DGX 24( 12) ,0GU24( 12) . 
*OGX J2( IS) , DGW32( lb) , DGX40C 20) . DGW40( 20) 

DATA DGX I 2/.9ai5b063424bT19, .90411T25637047S,. T6990 2674 l 9 4 30 5 , 
*.537317 95 42 36 617, . 36 7331 49 39 9.3 130, .125 23340 3511469/ 

DATA DGW 12/. 4717 S 1363365113C-1*. 10 6939325995313, 

*• 160073323543 346».20316742672306S».2334925J653B355» 

*.2491470493 1 340 3/ 

DATA 0GX16/. 9394 00934991650, .944575023073233, . 36 56 1 1 202 3873 32 * 
*.7554 04 40 3 355003, . 6 l 7 3 7b 2 4 4 4 0 2 64 4 ,. 4 580 1 6 77 76 5 7 2 2 7 , 

*.2Slb 03550779259, .950125 09 33763740-1/ 

DATA DGW 16/. 27 l 5 24 591, 11 754 ID- l*. 6 2 25 352 39 3 36 4790-1* 

*. 95 1 5 35 l lb 32492 3 D-l». I 2462 397 1255534 *. 149 59 59333 16 577, 

*. 1691 56 51939 5001*. 1326 0 341 504 4924, . 189450610455063/ 

DAT A DGX 24 /.9951 3 721 99-97021 , .974 723555971309, 

*.9 33274 5 5200 27 3 5*. 8 1641 552700 4401*. 3200 0193597 390 3* 

*. 740 124191573554*. 64 809 36 5 1936976 *. 545421471 38 38 40, 

*.4 3 37 >350762 6045, .31504 26 7*696163*. 191 11336 74736 16* 
*.64056892862b05bD-l/ 

DATA DGW24/. 12 34 1 2297 99937 2 C- 1, .285 3138 A623 9 3 370-1* 

*.44 27 74 3331 74 1 93 0-1 ,. 59 29 35 36 91 54 36 3D- 1, . 73 34 6 43 141 10 80 30-1 , 

*. 361 90 161 5 319531 D- l, .97 61365210411 3 9 D-i*. 107444270115 96 6, 

*. I 15^ 156b 305 372b, . 12167047 ’92 730 3. . 125837456346528, 

*. 12 79 33195 146 7 52/ 

DATA DGX 32/. 99 726 3 361 3494.3 2, . 985611511545268 , 

*.964762255537506 , .934906075937740*. 396321 155766052* 
*.8495670137325 TO, .794 48 3795 9b 7942, .7321321187402 90, 
*.6fa304426b9302l5,.5377l5757240762,.506899908932229, 

* .421 351276 l 306 35 ,. 3 11 36860 223 2123, 

*.239 2 37362 2 52137, .14447196153279b, .4330 7665 6877383 0-1/ 
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data rv.w 3 2/. 7C 10 b 1 coo 44 to l on- 2, . 162743947 3390570 - 1 , 

9 » 2 E 3 '‘.VS I. G 30 92 *> 2 l 0 - l * • 3 4 2 7 3 9 6 2 9 1 3 0 2 1 4 0- 1 * • 4 2 3 3 5 9 9 8 0 2 2 2 2 6 7 D - 1 * 

* .309 *9059262 J7-, ’0-1 , .93634 39 397 0 5 3 550-1 , .S33 22 2227 7b 36 180-1 , 

*. 7 2 3 9 >7 , . 78 193 995 7.0 7 070 3 0-1 9 3 3 1 1 92 9 2 2694 ft AD- 1 , 

3 .3 76 5 209 303443 3.30-1 ,. 'U l 7 307.9 645 76 3 9l)- l , .9 39 44 399 0 90804 60-1 » 

9 ",b 3 372007927 490-1 ,.9-, 94 00 9951 4 7 2790-1/ 

OaTA O-.K 4.0/. 99 02 17709 7105 59, . 990 724 2 386994 >7, 

S. 9 77, Ti *949 90 J774,. 95 7 4 16319213792, .93281290 827. 96 77, 

3 .9O2O ‘0.9O69S0I3 74 , . 9 >>5 95 9 50 3212260, . 9246 122 300 3 3 312, 

3. 77 fll. 15*. 514 2»51*,. 72731925?. 10 *9 2 7, .4719^0654614100, 

3 .61 254 39 090^7 9 00 , .5-, 9 4-. 7 1 2". O'* 5 l 2 3 ,. 4 9 30 7 59 3 1 6 3b l 7 9 , 

3 .4 1 37 79204 3 7 ) i ,. »4 19940 » "2 .7 5 9, . 26A 1521 3 500 72 54 , 

3. I *2~ 9 70 00 70 1 5 7 1 , . 1 140 0 40 7-J *»7 5 25 5 , . 30 7 7 24 l 7 00605 08 0- 1/ 

9 A I a DGW 40/. 4 5 2 12 7 ’.1 9 V. 3 3 1* 0-2, .1049820453115200-1, 

3 . 1 442 1000 35 l wr: 7 40- 1 ,. 22 24 5 S4 9 1*4 16 7 00- 1 ,. 979 37 006 90002 340-1 , 

3 . 3 344 019 52-32 -4 7.- 1 f>— 1 , . 3 3 792167 9744 7 20 0-1 ,.4 3 8 7 0 90 8 1 0 567 3 3 0 - l , 

3 .400949074 1507 2 2 0-1 ,. 5 322’ 44^93 39 36.00- 1 ,. 57 4 3 974 9099 39160- l , 

3.41 3042424929299.9- 1 ,. 64 904 ; l 345660100-1 , .679120459 1 52339D-1 , 
3.7C6116473912. 3030-1, .72 39 6 492395 80 41 0-1, .74723169057 96030-1, 

3 . 76 11. 3 34 l 90 062620-1 ,- 7 70 393 101642480 0-1 ,. 77505 94 7 97 042400-1/ 

e x t e r nal dfun 

C CALCULATE THF 12 POINT AND 16 POINT INTEGRAL ESTIMATES 
DINT I 2 = OGOUAO( OLOWE A , O'JPPE- , OF UN, 6, DGX 12, 0GW12) 

DINT 16 = OGOUAO( DLOUE0 , OUPPE* , DFUN, 8 , OGX 16, DGM16) 

C DETERMINE WHETHER CONVERGENCE WAS OBTAINED 

CALL CONVTI D INTI 2, DINT 16,09 FLER, 0A8SEH , DERR , ICONV) 

I F { ICONV .Nt . 1) GOTO 10 

C CONVERGENCE - SET UP VALUES, and return 
I F L A G = 1 
DINT=0INTt6 
RE TURN 
C 

C NO CONVERGENCE - CALCULATE THE 24 POINT INTEGRAL ESTIMATE 
10 DINT2 4 = OGOUAO( OLOWFR, DUPPE 'R, DFUN , 1 2 , DG X 24 , DGW 2 4 ) 

C DETERMINE WHETHER CONVERGENCE UAS OBTAINED 

CALL CONVT< 0 INT 16 ,01 NT 24 , ORELS* , DA USER , DERR , IC ONV) 

I F ( I C ON V .NE. 1) GOTO 20 . . 

C CONVERGE VC E - SET U a VALUES, AND RETURN 
I F L A G = 2 
niNT=0INT24 
RETURN 

c 

C NO CONVERGENCE - CALCULATE THE 32 POINT INTEGRAL ESTIMATE 
20 DINT 32 = DGOUAOI OLOWER , OUPPFR, OF UN, 16 ,DGX 32, 0GW32) 

C !) E T F A M I N c WHFTncs CONVERGENCE WAS OBTAINED 

CALL CONVT( DINT >4 ,0 INT 3 2, ORELER, PARSER , DERR , ICONV) 

IF( IC.1NV .NE . 1) GOTO 30 
C CONVERGENCE - SET UP VALUES, AND RETURN 
I F L A G - 3 
DINT rOINT 32 
return 
C 

C NO CONVERGENCE - CALCULATE THE 40 POINT INTEGRAL ESTIMATE 
30 01 NT 4 Ori'G JU AD( DLOUf R , DUPPER , D F UN , 20 , OG X 4 0 , 0GU4 0 ) 

C «E TFRMIN" WHETHER CONVERGFNCF WAS OBTAINED 

CALL CSINV T( 0 INT 32, DINr40,9RFi.ER,DABSER, DERR, ICONV) 

IF (ICONV .NE. 1) GO T 0 4 0 
C CONVERGENCE - SET 'JP VALUES, AND RETURN 

I c l a . ; - 4 

DINT =ri INT 40 I. 

RETURN 

C 

C NO CONVERGENCE 
40 ifla.; = o 

DINTrO INT 40 
R ? TIJPU 

END 

SU HR OUT INF CONVT ( DOLO, JMFW , DA EL E •* » D« 3 SE R , DE R » , IC ON V ) 

C SUBROUTINE TO TEST CONVERGENCE FDR THE SUBROUTINE 3GAUSS. DOLD ANO 
C ONEW ARE THF PREVIOUS a NO PR “SENT ESTIMATES OF TH C INTEGRAL, 
c 9-;r: mt paper narFn aug 22 . 

c 

IMPLICIT PE AL3.0 (0) 

rcoNv=i 

DERR=I)ARS ( ONE W- DOLD) 

C IF THE PRESENT INTEGRA*. ESTIMATE IS EOUAL TO ZERO, PROCEED TO THE 
C ABSOLUTE ERROR TEST (AFTER making SURE THAT 0A9SER .NE. 0). 

I F ( 0 N E W .fO. 0 • DO ) GOTO 10 

C ESTIMATE NUT CJUAL TO ZERO - CHECK THE RELATIVE ERROR CRITERION 
IF ( DEBR/flARSf ONEW) .LF. DRELER) RETURN 
C NO RELATIVE fhrcjr CONVERGENCE - IF DABSER = 0, THEN THERE IS NO 
C CONVERGENCE. 

10 IF ( PARSER . EO. 0.00) GOTO 20 

c check the absolute error criterion. 

IF (DERR .LE. 0 A 3 S L R ) RETURN 
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C Ml CONVE 'Sr NC‘ : 

<0 IC'IMVO 
R E TURN 
END 

0 0 U R L E PRECISION FUNCTION OGQU A D ( DL OWE a , OU P P £ R > DF UN , NP T 02 , OGX , 
o DGW) 

C FUNCTION SUBPROGRAM TO ESTIMATE THE INTEGRAL FROM DLOWER TO OUPPER 
C OF THE FUNCTION DFUN USING GAUSSIAN QUADRATURE WITH 2SNPT02 POINTS. 

C ARRAYS OGX ANO DGU CONTAIN THE NORMALIZED ARGUMENTS ON THE INTERVAL 

C (0,1) (UITH f)GX(l) ThE POINT CLOSEST TO 1) AND ThE RESPECTIVE 
C WEIGHTS. 

C REP; NUMERICAL METHODS, RY HORNRECX, SECTION fl.<* P. 15<». 

C MY PAPER DATED AUG22. 

C 

IMPLICIT REALS* (D) 

REAL** DGX(NPT02) ,0GW(NPT02> 

C CALCULATE The CONSTANTS USED IN FORMING the arguments of ofun. 

DC 1 = ( DUPPFR4-0L0WER)/ 2.1)0 
DC 2 = DC 1- DLOWfcR 

C NOW PERFORM THE INTEGRATION 
DGQU A D = 0 • DO 

DO 10 I-1,NPT02 „ ... 

10 DG QUA D = DGQU AI)f DGU ( I ) S ( 0 FUN ( DC 1 ♦ DC 2« DGX ( I ) )*DFUN( DC l -DC 2« DGX { I ) ) ) 
DGQU A D = DG Q'JA OS DC 2 . 

RETURN 

FND 
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FORTRAii SCURCC CODE FOR CALCULATIiiG 
FORCE AT TIME T=50 MICROSECONDS 


' rL r :: FORWU FORTRAN VI. 1 

PP.OG?<A" TC ESTIMATE ThF TOTAL FO = CE pET«' 7'I THE COIL AND TH C SKIN 
r-Y CALCULATING < 10**71/4 -IM£S THF f S r I M a T £ of ThF I N T tf G* A L OF 
BR*BR*R - BZ*3Z*R 
BETWEEN P=.C1 AN D R = 1.T INCHES. 

ref: my paper dated jan n •«<». 

REAL BPSONC 11) / .0 17014, . ?•? 6 44 <5 , . 45 ;, 5 52 , . 50 3 3 3*.,. 4 AIT 93, . 2 3 '•» 1 <» 3, 
P.tC/34 3 7,.04flIA0,.0 2SI20,.G'I44 7«,.0CE‘ , :'3/ 

REAL n Z SO N( 111/. 2244 64, .213477, .142* 2’, .05 27 IS, .0 33 232, -.' 31324, 

* — • 047 3 S 4 • - »0 17 1 9*»» — • 0Cb4 34 » — • C02**4 1 t “• 1629/ 

Pf "L DASONSt 11 > »CR50NS( 1 1) » DR50NS( 11 >» 3Z 50NS( 11),CZS0NS(11), 
*CZSCNS( 11) 

REAL AA{ll)/.ni,.2,.4*.6».?*l.,l«2»1.4»1.6,l.ft»2./ 

REAL RZ(II)/. 01, . I, .3».5.« 7, .4,1. 1,1. ’,1.5, 1.7, 1.9/ 

REAL MONTH/* JAN* /,C ATE/* 17 •/ 

DATA C0NST/2.SE6/, T0MTRS/.0254/, TOLB S/ .2248/ 

N=ll 

SQUARE THE CALCULATED B FIELDS 

00 1C 1=1,11 

3 A "50 N( I)=BA50N( I)*3P50N< I) 

1C BZ50NC t) = BZ50N( T)*BZ50N( I) 

CHANGE THE RADII AT WHICH T HP 8 FIELDS WERE CALCULATED FRO** INCHES 
TO METERS 

DO 15 1=1,11 
PH(I j=RR(I)*TOMTRS 
15 RZ( I ) = P Z ( IJRTOMTRS 

Eir SPLINE FUNCTIONS TO 8R*BR AND B Z*3Z 

CALL SPLINE (N, RR,BP50N,PR50NS ,CA SONS ,0 R50NS1 
CALL SPLINE(N»RZ,BZ50N,BZSCN5 ,CZ50NS ,0 7 SON S) 

Pcpfqrm THE INTEGRATIONS 

R I NT SO =0. 

Z IN TSO =0 . 

LOOP TO ADD THE CONTRIBUTIONS TO THE INTEGRAL BETWEEN KNOTS 
00 20 1=1,10 
TP 1 = 1* l 

PINTS0 = RINT5C*SPLINT(I,RP. (IP1 ),R.R( I) ,1 R GON( I ) , BASONS C I) ,CR50NS( I) , 

* DR«ONS(in 

>0 ZINTF0=2INT507SPLINT(t f RZdPl )»RZ( t ) ,3 Z50NC I ) , BZSONSC I) ,C7S0N5( I) , 
o DZSONS(I)) 

DONE EXCEPT FOR MULT I PL IC AT I ON BY CONSTANT 

FSDMET=CONST«(R INT50-Z IN TSO > 

Fr0LSS=F5CMET=70LBS 

UR I T E ( 6 , 1000 ) MONTH, DATE ,F50MET,FS0LBS 
1 '* 7 0 FORMAT (*10UTPUT FROM FORWU FORTRAN V 1. 1 * , 1 OX , A 4, A 4 ,// • TOTAL FORCE 

* AT SOUSEC IS * , F 9 . 4 , * NEWTONS*, /, T25, •= *,F9.4,* POUNOS*,//) 

STOP 

END 

FUNCTION SPLINT (I, PUP, RLO,F,E,C,D) 

FORTRAN SUBPROGRAM TO EVALUATE THF INTEGRAL OF F(R) o R 

OV fs the INTERVAL RLO .LE. P .LE. RUP WITH 

1 = SPLINE FUNCTION INDEX COP P E SPOND IMG TO THIS INTERVAL 

PUP = UPPER INTEGRATION LIMIT (A KNOT! 

RLO = LOW c A INTEGRATION LIMIT {A KNOT) _ 

F = F(=LO) 

P = 2(1) IN THE SPLINE INTERPOLATING "'LVNQMIAL OF ECO 

C = C( I j IN the SPLINE . .. 

C = CCt) IN THE SPLINE . .. 

A S >U *' * S A CUBIC S=LINE INTERPOLATING. POLYNOMIAL AS DESCRIBED Tv 
'•COMPUTER METHODS ‘■OP MATHEMATICAL COMPUTATIONS", BY FOASYThF C T. 

OL. 

hfe: »y pop er dated jan 14 ps. 
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■»Lr. : = a L(!®''LC 
RL r, .‘ ==Lc: j ^Ln 
p.LC^RLCT^RLQ 
PLO's’LtXfS'KLO 
RU?.’-iiUP* t, JP 

pupj=rup2*bup 

PUP4=EljP3*PUP 
RUPf = « UP4P l UP 
CALCULATE FIRST T E Rri 

S»LTMT = ( ( ( DPP.CO-C )o«LH4 B)PPL0-P »P (»LP! - TUP2)/2. 

act seconc tep* 

splint =SPLINT+( O .“DoPLQ-2. «C )*P.LCUP )" (PUP 3-RL03 )/3. 
Al'C THIRD TE°« 

SPLINT = SPLINT-{ 3.PDPRL 0-C )*(RUP4-RLO<») /<* . 

ATI) c 0U 3 TH TE1HJ DONE 

SPLINT=SPL INT*DP( RUP5-flL0=)/5. 

RETURN 

FND 


ORIGINAL 
°f Pool? 


PAGE IS 
QUALITY 


SUSP OUTINE SPLINEC N » Xt Y, B t C ,D ) 
CCPItO SUBROUTINE SPLINE FOLLOWS 
INTEGER N 

REAL X(N)'Y(N) f E(N)«C(N)»C(N) 
INTEGER NV1, I3t I 
PEAL T 
N'n = N-l 

IF(N.LT.2)RETURN 
IF(N.LT.3)GO TO SO 
D(1)=X(2)-X(1> 

C(2)=( Y(2)-Y<1) )/D(l> 

M 15 I = 2 » N *» I 

0(1) = x ( 1*1 )-X(I ) 

P( I) = 2«e(D(I-l)»D(I )) 

C ( r ♦ 1 ) = (Y ( 1*1 )-Y( I) )/D<I ) 

cci ) = cci*i)-cm 

10 CONTINUE 
B ( 1 ) =-D ( l ) 


l 1 


n^ic vtu ' i> 


C(N)=0. 

IF(N.E3.3) GO TC 15 

C(1)=C(J)/(X(4)-K(2)) - C(2)/(X(3)-X(i ) ) 

C(N)=C (N-l )/(X(N)-X(N-2) ) - C(N-2)/{X(N-l)~X(N-3>) 
C( 1 )=C Cl )<• 0(1)002/ <XC«*)-X(1)) 

C(N) = -C<N)~ D<N-I)«e2/(X(N)-X(N-?) ) 

DO ?C I = ?, N 

T=ofi-i>/B< r-i) 


ni 

Cn CONTINUE 

C(N)=C (N)/3(N) 
no <o IF.= 1, ■|“1 


I=N-T3 

c ( n=(ca )-cci >*c ci>i >)/e( i) 

•0 CONTINUE 

D(N> = ( Y(N)-Yf ft!Ml ))/l (N«1 ) ♦ rC'IVl )■>{ C(NWl) * 2.*C(N)) 

” ?!!( : * ,,!,Mr<, *»* 2 -' c “» 

C{ I ) = 3 »"C ( I ) 

40 CONTINUE 

C(N)=3.*C(N) 

DjN)sO(N-l) 

RETURN 

50 PC I ) - ( Y (2 )— Y ( 1 ) )/(X(2)-X(i)) 

C(l)=C. 

oil )=n . 

B(2)=E<1) 

C ( 2 ) *0 • 

D{ 2)=0 . 

RETURN 

END 



APPENDIX F 


NOTATION AND LIST OF SYMBOLS 

Unless stated othervi.se, the rationalized MKS system is used in 
all equations and calculations. Equations are numbered consecutively, 
beginning with (1) in each chapter. Only those equations that are 
referenced in the text are numbered. Figures are identified by tvo 
alphanumeric characters. The first character is a number, indicating 
the number of the chapter in which the figure appears. The second 
character is a letter, alphabetically identifying the order in which 
the figures appear. 

A chapter by chapter summary of the symbols used in this disser- 
tation, listed in the order of their appearance, follows. 

Chapter One 

Z(oj ) Terminal coil impedance (page 7) 
co Radian frequency (Fourier transform) variable (page 7) 

i(t) Time domain coil current (page 7) 

I (co) Fourier transform of i(t) (page 7) 

E^(z, r, t) Azimuthal real space component of the electric in- 
tensity (page 7) 

z Axial coordinate in cylindrical coordinate system 

(page 7) 

r Radial coordinate in cylindrical coordinate system 

(page 7) 
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H, (z, r, t) Radial real space component of the magnetic intensity 
< page 7 ) 

H, (z, r, t) Axial real space component of the magnetic intensity 
(page 7) 

f(t) Separation force magnitude between the coil and metal 
target (page 7) 

r 1 Impulse delivered to metal target by the coil (page 7) 


Chapter Three 

Ri Inner coil radius (page 12) 

R* Outer coil radius (page 12) 

h Coil axial width (page 12) 

d Metal target thickness (page 12) 

$ Azimuthal coordinate in cylindrical coordinate system 

(page 12) 

^ Permeability of free space (all materials considered in 

this dissertation are non-magnetic) (page 14) 

o - Conductivity of metal target (page 14) 

j Principal square root of -1 (page 15) 

J, (x) Bessel function of the first kind of order n (page 15) 

X Hankel transform variable (page 15) 

A 

V Total Fourier-Hankel space transmission line voltage 

(page 16) 

I" Total Fourier-Hankel space transmission line current 

(page 16) 

Fourier-Hankel space transmission line complex propaga- 
gation coefficient (page 17) 

♦ Fourier-Hankel space transmission line phasor current 

corresponding to propagation in the direction of in- 
creasing z (page 17) 

A 

I«- Fourier-Hankel space transmission line phasor current 

corresponding to propagation in the direction of de- 
creasing z (page 17) 
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V 0 . Fourier-Hankel space transmission line phasor voltage 

corresponding to propagation in the direction of in- 
creasing z (page 18) 

Vo- Fourier-Hankel space transmission line phasor voltage 

corresponding to propagation in the direction of de- 
creasing z (page 18) 

Z s Fourier-Hankel space transmission line characteristic 

impedance (page 19) 

I Fourier space phasor current (page 19) 

g Distance between metal target and the current sheet 

closest to target (page 20) 


z 0 Axial coordinate of current sheet closest to target 

(page 20) 

Zi Axial coordinate of center current sheet (page 20) 

z* Axial coordinate of current sheet furthest from 

target (page 20) 


Chapter Four 

V Fourier space phasor voltage (page 23) 

Z« Fourier-Hankel space air transmission line characteris- 

tic impedance (page 28) 

X Fourier-Hankel space air transmission line complex propa 

gation coefficient (page 26) 

Z, Fourier-Hankel space metal transmission line character- 

istic impedance (page 26) 

Fourier-Hankel space metal transmission line complex 
. propagation coefficient (page 26) 

Z(0) Fourier-Hankel space impedance seen looking into metal 
transmission line at z=0 in Figure 3C (page 26) 

/O^ Fourier-Hankel space current reflection coefficient 

(page 26) 

p 6 Fourier-Hankel space voltage reflection coefficient 

(page 28) 

Fourier-Hankel space total voltage at coil side face of 
metal target (page 30) 


V, 
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V, Fourier-Hankel space total voltage on face of target 

opposite coil (page 30) 

I, Fourier-Hankel space total current at coil side face of 

metal target (page 30 > 

If Fourier-Hankel space total current on face of target 

opposite coil (page 30) 

F Vector force on metal target (page 32) 

£ Permittivity of free space (all materials considered in 

this dissertation have a relative permittivity of 1) 
(page 32) 


Chapter Five 


r» Radius of inner loop of wire in field measuring plate 

(page 38) 

r« Radius of outer loop of wire in field measuring plate 

(page 38) 

h Distance between two loops of wire of the same radius on 

front and back sides of field measuring plate (page 38) 


Chapter Six 


R(oj) 
L( to) 
V. 

SM 

V„ (to ) 


C 

Re c 

\ } 

Io 


Real part of the total coil impedance (page 47) 
Inductance of coil (page 47) 

Initial voltage on energy storage capacitor prior to 
discharge (page 47) 

Dirac delta distribution (page 47) 

Fourier transform of energy storage capacitor voltage 
(page 47) 

Capacitance of energy storage capacitor (page 48) 

AC winding resistance of the coil (page 48) 

Symbol denoting *real part of ■£ "]- ■ (page 48) 

Instantaneous current in coil when clamp diode across 
energy storage capacitor begins conducting (page 51) 
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Appendix C 

H«(x) Struve function of order n (page 80) 

T» (x) Chebyshev polynomial of the first kind of order n 
(page SO) 

Y»(x) Bessel function of the second kind of order n (page fll> 

Partial sum associated with the Jacobi series expansion 
of H„ (x)-Y„ (x) (page 01) 


A„ (x! 
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